Numerical Scheme for Solving the 2D-Temperature Field

Boundary conditions

Boundary conditions both for conduction and convection are needed. For convection boundary conditions can be in the form of a heat transfer number and an ambient temperature on all the outer boundaries. The way the contributions from these boundary conditions are inserted into the equations is very similar to how the internal conductances (see the chapter "The dynamic 2D temperature integration to read more about the internal conductances) are computed and inserted into the equations.

On the outer diameter of the cylinder the computation is closed with thermal boundary conditions. Enough rock must be included in the computation so that these boundary conditions do not influence the temperature computation in the well to any large extent. In case the simulation is done with a riser, the cylinder is extended with the sea water surrounding the riser. The outer part of the cylinder, sea water, is treated as a dummy area ('Void material'), and the given thermal boundary conditions are applied directly to the riser surface.

When having a riser and a wellhead, they will be surrounded by water and air. No simulation will be done in the water and air material ('Void material'). Grid cells with the void material are treated in a special way in the program:

• The cell temperature is simply set to zero ($T = 0\degree\text{C}$)

• The boundary conditions given on the boundary of the cylinder are translated into the boundary of a Void/NoVoid material (Riser surface/wellhead surface as shown in the figure). The same also applies to the top boundary conditions that are translated down to the sea bottom/top of the wellhead.

When simulating a production case, the temperature of the produced fluid at the depth where it is produced is considered as a constant temperature that is the same as the geothermal temperature. When simulating an injection case, the injected fluid temperature can be chosen by the user. In circulation or drilling mode, when a mudpit is used, it is possible to set up the mudpit temperature.

Initialization

Each gridcell needs to be initialized with two important state variables:

• The central temperature in the cell. True vertical depth as a function of measured depth and initial temperature as a function of true vertical depth are input to the program as tabular functions. The initial temperature can therefore be computed for each grid cell.

• The material in the cell. This material evolves dynamically when successive temperature simulations are done during the construction of the well.

Time Integration

After assembling all the contributions to the energy balance as detailed above, a set of ordinary differential equation is obtained. They can be written in matrix form as:

$$D \frac{\partial{\vec{T}}}{\partial{t}} = B \vec{T} + \vec{b}$$

$D$is a diagonal matrix (the accumulation term) and $B$ is a band structure with five bands, the diagonal and four bands for the four neighbors. By integrating this equation system between the time $t_0$ and the time $t_0 + \delta_t$, one gets:

$$D \frac{\vec{T}-\vec{T_0}}{\, \delta_t} = B {\vec{T}}_{\text{mean}} + \vec{b}$$

$T_{\text{mean}}$is the time averaged temperature in the time interval. By approximating this temperature with the known starting temperatures $T_0$, we get the so-called explicit Euler integration scheme. By approximating the temperature with the end point temperature $T$, we get the so-called implicit Euler scheme. The explicit scheme has the advantage that the solution can be obtained just by matrix multiplication, but it has limited stability properties. The implicit scheme is unconditionally stable, but requires us to solve a set of equations. The implicit Euler scheme is used to integrate the equations. Rearranging terms in the equation system, one gets:

$$\left(\frac{D}{\, \delta_t} - B \right) \vec{T} = \frac{D}{\, \delta_t} \vec{{T}_{0}} + \vec{b}$$

The matrix is positive definite, but not symmetric due to the convective terms. A standard direct solver for band matrices is used to solve the equation system, providing us the new temperatures at the end of the time step.