Skip to main content

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βˆ’zr-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:

d dt∭u(T,p)ρ(T,p) dV⏟AccumulationΒ termΒ =∬qn dA⏟ConductiveΒ termΒ +∬en dA⏟ConvectiveΒ term+∭sn dV⏟SourceΒ term\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}



  • ρ\rho is the mass density
  • uu is the specific internal energy (per mass)
  • qq is the conductive energy flow
  • ee is the convective energy flow
  • ss 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, hh, as:

ddt∭u(T,p)ρ(T,p) dV=ddt∭(h(T,p)βˆ’pρ)ρ(T,p) dV\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

TT and pp 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 pp is also assumed. Then this accumulation term can be written as:

dhdtρVβˆ’dpdtV=dhdTρVdTdt=cpρVdTdt\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 are computed in the 1D flow computation:

(βˆ‘iΞ±icpiρi)VdTdt\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: qβƒ—=βˆ’kβˆ‡T\vec{q} = -k\nabla T

where the conductive heat flow is proportional to the temperature gradient. The proportionality factor, kk , 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, QuQ_u , QbQ_b , QlQ_l , and QrQ_r can be expressed as:

Qr=Cr(Trβˆ’Tc)Q_r = C_r\left(T_r-T_c\right)

Ql=Cl(Tlβˆ’Tc)Q_l = C_l\left(T_l-T_c\right)

Qu=Cu(Tuβˆ’Tc)Q_u = C_u\left(T_u-T_c\right)

Qb=Cb(Tbβˆ’Tc)Q_b = C_b\left(T_b-T_c\right)

where CrC_r, ClC_l, CuC_u, and CbC_b are the conductances between the grid cells in the different directions. TrT_r, TlT_l, TuT_u, and TbT_b are the temperatures fo the neighboring cells and TcT_c is the temperature of the central cell.

When doing this, it must be taken into account that in general the thermal conductivity kk 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 CrC_r, ClC_l, CuC_u, and CbC_b 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 zz direction are given by the thermal conductivity kk, the area AA where the conductive flow goes through and the distance between the central point and the boundary of the cell:

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

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

Cri=k2Ο€dzln⁑(r2r1){C}_{r}^{i} = k \frac{2 \pi dz}{\ln \left(\frac{{{r}_{2}}}{{r}_{1}} \right)}

where r1r_1 is the average radius of the cell and r2r_2 is the radius of the external boundary of the cell. A similar expression is obtained for CliC_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:

Cr=CriCli+1Cri+Cli+1C_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:

Qr=Cr(Trβˆ’Tc)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.

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:

Ξ”U=Fu δt\Delta U = F u \, \delta_t

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

Ξ”V=F δtρ\Delta V = \frac{F \, \delta_t}{\rho}

Adding the internal energy transferred and the pΒ dVp\ dV work done by the upstream volume on the downstream volume, it is obtained:

Ξ”U+pΞ”V=Fu δt+pF δtρ=F δt(u+pρ)=F δth\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Β dVp\ dV 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Β dVp\ dV 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 (h(T,p)βˆ’gd+0.5v2)F\, \left(h \left(T,p \right) - g d +0.5 v^2 \right)

As dd is the True Vertical Depth, the equation is written with a minus sign. TT and pp should really be taken at the boundary. vv is the fluid velocity. For reasons of numerical stability, an upstream approximation is used, with the upstream values of TT and pp:

F (h(Tup,pup)βˆ’gd+0.5v2)F\, \left(h \left({T}_{\text{up}} , {p}_{\text{up}} \right) - g d +0.5 v^2 \right)

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

F(h(Tup,pup)βˆ’gd+0.5v2)=F(h(Tup0,pup)+βˆ‚hβˆ‚T(Tupβˆ’Tup0)βˆ’gd+0.5v2)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 TupT_{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 (hmixh^{mix}, cpmixc_p^{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.