Motivation

Assume the well is producing  of water,  of oil and  of gas as measured daily at separator with pressure  and temperature . 

While moving down to depth  along the hole the wellbore pressure  will be growing due to gravity of fluid column and friction losses emerging from fluid contact with inner pipe walls . 

Wellbore temperature  will be also varying due to heat exchange with surrounding rocks.

The volume shares , occupied by different phases will be varying along hole due to along-hole pressure-temperature variation, phase segregation and phase slippage.

Definition

Mathematical model of Multiphase Wellbore Flow predicts the temperature, pressure and flow speed distribution along the wellbore trajectory with account for:

• tubing head pressure which is controled by gathering system or injection pump
  
  
• wellbore design (pipe diameters, pipe materials and inter-pipe annular fillings)
  
  
• fluid friction with tubing /casing walls
  
  
• interfacial phase slippage
  
  
• heat exchange between wellbore fluid and surrounding rocks via complex well design
  
  
  

Consider a 3-phase water-oil-gas flow: .

The -phase volumetric flow fraction ( also called phase cut or  input hold-up or no-slip hold-up ) is defined as:

\gamma_\alpha = \frac{q_\alpha}{q_t}

where  – volumetric flow rate of -phase and  is the total volumetric fluid production rate:

q_t = \sum_\alpha q_\alpha = q_w + q_o + g_g

In multiphase wellbore flow each phase occupies its own area  of the total cross-sectional area  of the lifting pipe. 

This area can be connected into a single piece of cross-sectional area (like in case of slug or annular flow) or dispersed into a number of connected spots (like in case of bubbly flow).

A share of total pipe cross-section area occupied by moving -phase is called an -phase in-situ hold-up and defined as: 

s_\alpha = \frac{A_\alpha}{A}

so that a sum of all in-situ hold-ups is subject to natural constraint:

\sum_\alpha s_\alpha = s_w + s_o + s_g = 1

When word hold-up is used alone it usually means in-situ hold-up and should not be confused with input hold-up or no-slip hold-up  which should be better called volumetric flow fraction.

The actual average cross-sectional velocity of moving -phase is called in-situ velocity and defined as:

u_\alpha = \frac{q_\alpha}{A_\alpha}

where  is the volumetric -phase flowrate through cross-sectional area .

The superficial velocity of -phase is defined as the 

u_{s \alpha} = \frac{q_\alpha}{A}= s_\alpha \cdot u_\alpha

The multiphase mixture velocity  is defined as total flow volume normalized by the total cross-sectional area:

u_m = \frac{1}{A} \sum_\alpha q_\alpha = \sum_\alpha u_{s \alpha} = \sum_\alpha s_\alpha \cdot u_\alpha

The difference between velocities of -phase and -phase is called interfacial phase phase slippage:

u_{\alpha_1 \alpha_2} = u_{\alpha_1} - u_{\alpha_2}

The multiphase fluid density  is defined by exact formula:

\rho_m =  \sum_\alpha  s_\alpha \rho_\alpha

where  – density of -phase.

The two-phase gas-liquid model is defined in the following terms:

u_m = u_{s g} + u_{s l} = s_g u_g + (1-s_g) u_l

The two-phase oil-water model is defined in the following terms:

u_m = u_{s o} + u_{s w} = s_o u_o + (1-s_o-s_g) u_w

The 3-phase water-oil-gas model is usually built as a superposition of gas-liquid model and then oil-water model:

Input & Output

Application

Mathematical Model

The multiphase wellbore flow in hydrodynamic and thermodynamic equilibrium is defined by the following set of 1D equations: 

\frac{\partial (\rho_m A)}{\partial t} + \frac{\partial}{\partial l} \bigg( A \, \sum_\alpha \rho_\alpha \, u_\alpha \bigg) = 0

\sum_\alpha \rho_\alpha \bigg[ \frac{\partial u_\alpha}{\partial t} + u_\alpha \frac{\partial u_\alpha}{\partial l}  - \nu_\alpha \Delta u_\alpha\bigg]  =  - \frac{dp}{dl} + \rho_m \, g \, \sin \theta - \frac{ f_m \, \rho_m \, u_m^2 \, }{2 d}

(\rho \,c_p)_m \frac{\partial T}{\partial t} 
 
-  \sum_\alpha \rho_\alpha \ c_{p \alpha} \ \eta_{s \alpha} \ \frac{\partial P_\alpha}{\partial t}  
 
+ \bigg( \sum_\alpha \rho_\alpha \ c_{p \alpha} \ u_\alpha \bigg) \frac{\partial T}{\partial l}
 \  =   \   \frac{1}{A}  \ \sum_\alpha \rho_\alpha \ c_{p \alpha} T_\alpha \frac{\partial  q_\alpha}{\partial l} 

where

Equations  –  define a closed set of 3 scalar equations on 3 unknowns: pressure , temperature  and mixture-average fluid velocity  .

The model is set in 1D-model with axis aligned with well trajectory :

The disambiguation of the properties in the above equation is brought in The list of dynamic flow properties and model parameters.

Equation   defines the continuity of the fluid components flow or equivalently represent the mass conservation of each mass component  during its transportation along wellbore. 

Equation  defines the motion dynamics of each phase (called Navier–Stokes equation), represented as linear correlation between phase flow speed   and pressure profile of mutliphase fluid .

The term  represents heat convection defined by the wellbore mass flow. 

The term  represents the heating/cooling effect of the fast adiabatic pressure change.

This usually takes effect in the wellbore during the first minutes or hours after changing the well flow regime (as a consequence of choke/pump operation). 

The term  defines mass-specific heat capacity of the multiphase mixture and defined by exact formula:

(\rho \,c_p)_m = \sum_\alpha \rho_\alpha c_\alpha s_\alpha

The in-situ velocities  are usually expressed via the macroscopic flow velocity  using the 

(\rho \,c_{pt})_p \frac{\partial T}{\partial t} 
 
-  \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \eta_{s \alpha} \ \frac{\partial P_\alpha}{\partial t}  
 
+  \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ u_\alpha \frac{\partial T}{\partial l}
 \  =   \   \frac{\delta E_H}{ \delta V \delta t}

Stationary Flow Model

Stationary wellbore flow is defined as the flow with constant pressure and temperature:   and  .

This happens during the long-term (usually hours & days & weeks) production/injection or long-term (usually hours & days & weeks)  shut-in.

The temperature dynamic equation  is going to be:

\sum_\alpha \rho_\alpha \ c_{p \alpha} \ u_\alpha \frac{\partial T}{\partial l}
 \  =   \   \frac{1}{A}  \ \sum_\alpha \rho_\alpha \ c_{p \alpha} T_\alpha \frac{\partial  q_\alpha}{\partial l} 

The phase temperature  is the temperature of the -phase flowing from reservoir into wellbore.

It carries the original reservoir temperature with heating/cooling effect from reservoir-flow throttling and well-reservoir contact throttling:

T_\alpha = T_r + \epsilon_\alpha \, \delta P = T_r + \epsilon_\alpha \, (P_e - P_{wf})

The discrete computational scheme for  will be:

\bigg( \sum_\alpha \rho_\alpha^{k-1} \ c_{p \alpha}^{k-1} \ q_\alpha^{k-1} \bigg) T^{k-1} - \bigg( \sum_\alpha \rho_\alpha^k \ c_{p \alpha}^k \ q_\alpha^k \bigg) T^k
 =   \sum_\alpha \rho_\alpha^k \ c_{p \alpha}^k \ (q_\alpha^{k-1} - q_\alpha^k) \, (T_r^k + \epsilon_\alpha^k \delta p^k )

where  is drawdown,  – formation pressure in -th grid layer,  – bottom-hole pressure across -th grid layer,  – remote reservoir temperature of  -th grid layer.

The -axis is pointing downward along hole with -th grid layer sitting above the -th grid layer.

If the flowrate is not vanishing during the stationary lift () then   can be calculated iteratively from previous values of the wellbore temperature  as:

(\rho \,c_{pt})_p \frac{\partial T}{\partial t} 
 
- \ \phi \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \eta_{s \alpha} \ \frac{\partial P_\alpha}{\partial t}  
 
+ \bigg( \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \epsilon_\alpha \ \mathbf{u}_\alpha \bigg)  \nabla P
 
+ \bigg( \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \mathbf{u}_\alpha \bigg) \  \nabla T 
 
 - \nabla (\lambda_t \nabla T) =  \frac{\delta E_H}{ \delta V \delta t}

u_\alpha = \frac{q_\alpha}{\pi r_f^2}

\bigg( \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \mathbf{u}_\alpha \bigg) \  \nabla T 
 =  \frac{\delta E_H}{ \delta V \delta t}

References

The list of dynamic flow properties and model parameters