Page tree

Synonym: = Volatile/Black Oil Reservoir Flow @model = Muskat - Leverett equation




Definition



Mathematical model of Volatile Oil reservoir flow predicts the temperature, pressure and flow speed distribution in reservoir with account for:

  • available historical data on surface flowrates and/or bottom hole pressure

  • available 3D geological model 

  • PVT and SCAL model

  • specific wellbore designs

  • gravitational forces

  • heat propagation

  • adiabatic and Jole-Thomson heat effects 


The Black Oil flow is specific type of the Volatile Oil flow with  R_v=0.


Mathematical Model



The Volatile Oil flow dynamics is defined by the following set of 3D equations:


(1) \partial_t \bigg [ \phi \ \bigg ( \frac{s_w}{B_w} \bigg ) \bigg ] + \nabla \bigg ( \frac{1}{B_w} \ \mathbf{u}_w \bigg ) = q_W (\mathbf{r})
(2) \partial_t \bigg [ \phi \ \bigg ( \frac{s_o}{B_o} + \frac{R_v \ s_g} {B_g} \bigg ) \bigg ] + \nabla \bigg ( \frac{1}{B_o} \ \mathbf{u}_o + \frac{R_v}{B_g} \ \mathbf{u}_g \bigg ) = q_O(\mathbf{r})
(3) \partial_t \bigg [ \phi \ \bigg ( \frac{s_g}{B_g} + \frac{R_s \ s_o} {B_o} \bigg ) \bigg ] + \nabla \bigg ( \frac{1}{B_g} \ \mathbf{u}_g + \frac{R_s}{B_o} \ \mathbf{u}_o \bigg ) = q_G (\mathbf{r})
(4) \mathbf{u}_w = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ ( \nabla P_w - \rho_w \ \mathbf{g} )
(5) \mathbf{u}_o = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ ( \nabla P_o - \rho_o \ \mathbf{g} )
(6) \mathbf{u}_g = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g - \rho_g \ \mathbf{g} )
(7) P_o - P_w = P_{cow}(s_w)
(8) P_o - P_g = P_{cog}(s_g)
(9) s_w + s_o + s_g = 1



(10) (\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) = \rho_{\rm inj} \, c_{\rm inj} \, T_{\rm inj} \, q_{\rm inj}({\bf r})\, \delta({\bf r})



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


The right sides of equations  (1) –  (3) suggest no sources of flow except the contacts between wells and reservoir which is specified by well models as boundary conditions (see below).






The Volatile Oil flow model simulates 3-component fluid : water, liquid hydrocarbon (called "oil") and gaseous hydrocarbons ( called "gas") that flow in 3 possible phases (water, gasified oil and free gas) and defined by the following set of equations:

(11) \partial_t \bigg [ \phi \ \rho_W \bigg ] + \nabla \bigg ( \rho_{Ww} \ \mathbf{u}_w \bigg ) = q_{mW}(\mathbf{r})
(12) \partial_t \bigg [ \phi \ \rho_O \bigg ] + \nabla \bigg ( \rho_{Oo} \ \mathbf{u}_o + \rho_{Og} \ \mathbf{u}_g \bigg ) = q_{mO}(\mathbf{r})
(13) \partial_t \bigg [ \phi \ \rho_G \bigg ] + \nabla \bigg ( \rho_{Go} \ \mathbf{u}_o + \rho_{Gg} \ \mathbf{u}_g \bigg ) = q_{mG}(\mathbf{r})
(14) \mathbf{u}_w = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ ( \nabla P_w - \rho_w \mathbf{g} )
(15) \mathbf{u}_o = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ ( \nabla P_o - \rho_o \mathbf{g} )
(16) \mathbf{u}_g = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g - \rho_g \mathbf{g} )
(17) P_o - P_w = P_{cow}(s_w)
(18) P_o - P_g = P_{cog}(s_g)
(19) s_w + s_o + s_g = 1




(20) (\rho \,c_{pt})_m \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}



Equations  (11) –  (13) define the continuity of the fluid components flow or equivalently represent the mass conservation of each mass component  \{ m_W, \ m_O, \ m_G \} during its transportation in space. 

Equations  (14) –  (16) define the motion dynamics of each phase, represented as linear correlation between phase flow speed   \bar u_\alpha and partial pressure gradient of this phase  \bar \nabla P_\alpha (which is also called Darcy flow  with account of the gravity and relative permeability).


Equations  (17) –  (18) define the hydrodynamic inter-facial balance between the phases with account of capillary pressure in porous formation  P_{cow}, \ P_{cog}. The key assumption is that capillary pressure at oil-water boundary is a function of  water saturation alone  P_{cow} = P_{cow}(s_w) and capillary pressure at oil-gas boundary is a function of  gas saturation alone  P_{cog} = P_{cog}(s_g)

In the absence of capillary pressure the inter-facial equilibrium simplifies and implies that all phases are at the same pressure at all times.  


Equations  (19)  implies that porous space is fully occupied by fluid at all times  \{ s_w, s_o, s_g \}.


Equation  (10)  defines the heat flow continuity or equivalently represents heat conservation due to heat conduction and convection with account for adiabatic and Joule–Thomson throttling effect.

The term  \frac{\delta E_H}{ \delta V \delta t}  defines the speed of change of  heat energy  E_H volumetric density.

In impermeable rocks ( \phi =0, \; \bar u_\alpha = 0) heat flow is defined by heat conduction only:

(21) \rho_r \, c_{pr} \frac{\partial T}{\partial t} - \nabla (\lambda_t \nabla T) = \frac{\delta E_H}{ \delta V \delta t}

The effective specific heat capacity of formation with multiphase flow is a simple sum of its components:

(22) (\rho \,c_{pt})_p  = (1-\phi) \rho_r \, \ c_{pr} + \phi \ (s_w \rho_w \, c_{pw} + s_o \rho_o \, c_{po} + s_g \rho_g \, c_{pg} )

The effective thermal conductivity of formation with multiphase flow is assumed to be a sum of its components:

(23) \lambda_{t} = (1-\phi) \ \lambda_r + \phi \ (s_w \lambda_w + s_o \lambda_o + s_g \lambda_g )

The term  \bigg( \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \mathbf{u}_\alpha \bigg) \ \bar \nabla T  represents heat convection defined by the mass flow. 

The term  \bigg( \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \epsilon_\alpha \ \mathbf{u}_\alpha \bigg) \bar \nabla P represents the heating/cooling effect of the multiphase flow through the porous media. This effect is the most significant with light oil and gases.


The term  \ \phi \sum_{a = \{w,o,g \}} \rho_\alpha \ c_{p \alpha} \ \eta_{s \alpha} \ \frac{\partial P_\alpha}{\partial t} represents the heating/cooling effect of the fast adiabatic pressure change. This usually takes effect in and around the wellbore during the first minutes or hours after changing the well flow regime (as a consequence of choke/pump operation). This effect is absent in stationary flow and negligible during the quasi-stationary flow and usually not modeled in conventional monthly-based flow simulations. 


The set  (11) –  (10) represent the system of 16 scalar equations on 16 unknowns: 

\{ T, \ P_w, \ P_o, \ P_g, \ s_w, \ s_o, \ s_g, \ u_w^x, \ u_w^y, \ u_w^z, \ u_o^x, \ u_o^y, \ u_o^z, \ u_g^x, \ u_g^y, \ u_g^z \},

which are all functions of time and space coordinates  (t, \mathbf{r}) = (t,x,y,z).


Expressing the molar densities with mass shares and phase density (see also "Volatile Oil Model") one gets:


(24) \partial_t \bigg [ \phi \ \rho_W \bigg ] + \nabla \bigg ( \rho_w \ \mathbf{u}_w \bigg ) = q_{mW}(\mathbf{r})
(25) \partial_t \bigg [ \phi \ \rho_O \bigg ] + \nabla \bigg ( {\tilde m}_{Oo} \ \rho_o \ \mathbf{u}_o + {\tilde m}_{Og} \ \rho_{g} \ \mathbf{u}_g \bigg ) = q_{mO}(\mathbf{r})
(26) \partial_t \bigg [ \phi \ \rho_G \bigg ] + \nabla \bigg ( {\tilde m}_{Go} \ \rho_{o} \ \mathbf{u}_o + {\tilde m}_{Gg} \ \rho_g \ \mathbf{u}_g \bigg ) = q_{mG}(\mathbf{r})
(27) \mathbf{u}_w = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ ( \nabla P_w - \rho_w \mathbf{g} )
(28) \mathbf{u}_o = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ ( \nabla P_o - \rho_o \mathbf{g} )
(29) \mathbf{u}_g = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g - \rho_g \mathbf{g} )
(30) P_o - P_w = P_{cow}(s_w)
(31) P_o - P_g = P_{cog}(s_g)
(32) s_w + s_o + s_g = 1


Substituting the values of mass densities and mass shares of fluid components (см. "Volatile Oil Model") and dividing each equation by density of corresponding component in standard conditions one gets the most popular form of Volatile Oil flow equations:





Initial Conditions



Initial temperature distribution is set as input:

T(0, \mathbf{r}) = T_0(\mathbf{r})


In case the simulation is performed over the undisturbed reservoir then initial temperature distribution is geothermal.


The initial condition on phase pressure, phase velocities and phase saturations is set by one of the following options: Equilibrium Start and Non-equilibrium Start.

Condition I – Equilibrium Start


Equilibrium Start means that flow was not happening before the start:  \{ \mathbf{u}_w = 0, \ \mathbf{u}_o = 0, \ \mathbf{u}_g =0 \} and correspondingly phase pressure \{ P_w, \ P_o, \ P_g \} and phase saturations  \{ s_w, \ s_o, \ s_g, \} were in stationary (not varying in time) conditions:

(33) \nabla \cdot \bigg ( \frac{1}{B_w} \ \mathbf{u}_w \bigg )_{t=0} = 0
(34) \nabla \cdot \bigg ( \frac{1}{B_o} \ \mathbf{u}_o + \frac{R_v}{B_g} \ \mathbf{u}_g \bigg )_{t=0} = 0
(35) \nabla \cdot \bigg ( \frac{1}{B_g} \ \mathbf{u}_g + \frac{R_s}{B_o} \ \mathbf{u}_o \bigg )_{t=0} = 0
(36) \mathbf{u}_w(0, \mathbf{r}) = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ (\nabla P_w(0, \mathbf{r}) - \rho_w \ \mathbf{g} ) = 0
(37) \mathbf{u}_o(0, \mathbf{r}) = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ ( \nabla P_o(0, \mathbf{r}) - \rho_o \ \mathbf{g} ) = 0
(38) \mathbf{u}_g(0, \mathbf{r}) = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g(0, \mathbf{r}) - \rho_g \ \mathbf{g} ) = 0
(39) P_o(0, \mathbf{r}) - P_w(0, \mathbf{r}) = P_{cow}(s_w)
(40) P_o(0, \mathbf{r}) - P_g(0, \mathbf{r}) = P_{cog}(s_g)
(41) s_w + s_o + s_g = 1


Condition II – Non-equilibrium Start


Non-equilibrium Start means that flow happening before the start:  \mathbf{u}_w^2 + \mathbf{u}_o^2 + \mathbf{u}_g^2 > 0  and correspondingly phase pressure  \{ P_w, \ P_o, \ P_g \} and phase saturations  \{ s_w, \ s_o, \ s_g, \} were in not in equilibrium:

(42) s_w(0, \mathbf{r}) + s_o(0, \mathbf{r}) + s_g(0, \mathbf{r}) = 1

pressure distribution  \{ P_w, \ P_o, \ P_g \} could be arbitrary providing the capillary constraints:

(43) P_o(0, \mathbf{r}) - P_w(0, \mathbf{r}) = P_{cow}(s_w)
(44) P_o(0, \mathbf{r}) - P_g(0, \mathbf{r}) = P_{cog}(s_g).


The phase velocities are initialized as:

(45) \mathbf{u}_w(0, \mathbf{r}) = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ ( \nabla P_w(0, \mathbf{r}) - \rho_w \ \mathbf{g} )
(46) \mathbf{u}_o(0, \mathbf{r}) = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ ( \nabla P_o(0, \mathbf{r}) - \rho_o \ \mathbf{g} )
(47) \mathbf{u}_g(0, \mathbf{r}) = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g(0, \mathbf{r}) - \rho_g \ \mathbf{g} )


In practice, the non-equilibrium conditions before the start  is usually a result of previous flow simulations for the same reservoir, sometimes using a different grid-structure.


External Boundary Condition 



The external boundary condition for the temperature is usually set by one if the two options:

External Temperature Boundary Condition I – Fixed Temperature

T(t, \mathbf{r}) |_{\Gamma_e} = T_e( \mathbf{r})

External Temperature Boundary Condition II – Fixed Heat Exchange

(48) \big( \mathbf{n}, \nabla T(t, \mathbf{r} \big) \big |_{\Gamma_e} = \zeta \cdot \big( T(t, \mathbf{r}) - T_e( \mathbf{r}) \big)

where  \zeta  – heat exchange coefficient at model boundary.


The external boundary condition for phase pressure, phase velocities and phase saturations is set by one of the two popular options:

External Pressure Boundary Condition I – Non-permeable boundary  \Gamma_e

(49) \big( \mathbf{n}, \ (\nabla P_\alpha(t, \mathbf{r}) - \rho_\alpha \mathbf{r}) \big) \big|_{\Gamma_e} = 0

where   \mathbf{n} – normal vector to the boundary  \Gamma_e and  \alpha = \{ w, o, g \}.

External Pressure Boundary Condition II – constant-pressure boundary  \Gamma_e

(50) P_\alpha(t, \mathbf{r}) \big|_{\Gamma_e} = P_i = const

where   \alpha = \{ w, o, g \}.

Well model



Well flow model (don't get confused with Wellbore Flow Model) simulates the flow at the contact between well and reservoir thus relating the sandface flow rates and pressure distribution in reservoir around the well:

(51) P_w(t, \mathbf{r}) \big|_{\Gamma_{WRC}} = P_o(t, \mathbf{r}) \big|_{\Gamma_{WRC}} = P_g(t, \mathbf{r}) \big|_{\Gamma_{WRC}}= P_{wf}(t) \big|_{\Gamma_{WRC}}

where  \Gamma_{WRC} is a well-reservoir.

Bottom-hole pressure  P_{wf}(t) \big|_{\Gamma_{WFC}} = P_{wf}(t,h)  at the contact (at depth  h along-hole) is set by one of the three popular conditions (traditionally called "Controls"):

  • Well Condition I – Pressure Control

  • Well Condition II – Liquid Control

  • Well Condition III – Oil Control


The list of dynamic flow properties and model parameters



(t,x,y,z)

time and space corrdinates ,

z -axis is orientated towards the Earth centre,

(x,y) define transversal plane to the z -axis

\mathbf{r} = (x, \ y, \ z)

position vector at which the flow equations are set

q_{mW} = \frac{d m_W}{dt}

speed of water-component mass change in wellbore draining points

q_{mO} = \frac{d m_O}{dt}

speed of oil-component mass change in wellbore draining points

q_{mG} = \frac{d m_G}{dt}

speed of gas-component mass change in wellbore draining points

q_W = \frac{1}{\rho_W^{\LARGE \circ}} \frac{d m_W}{dt} = \frac{d V_{Ww}^{\LARGE \circ}}{dt} = \frac{1}{B_w} q_w

volumetric water-component flow rate in wellbore draining points recalculated to standard surface conditions

q_O = \frac{1}{\rho_O^{\LARGE \circ}} \frac{d m_O}{dt} = \frac{d V_{Oo}^{\LARGE \circ}}{dt} + \frac{d V_{Og}^{\LARGE \circ}}{dt} = \frac{1}{B_o} q_o + \frac{R_v}{B_g} q_g

volumetric oil-component flow rate in wellbore draining points recalculated to standard surface conditions

q_G = \frac{1}{\rho_G^{\LARGE \circ}} \frac{d m_G}{dt} = \frac{d V_{Gg}^{\LARGE \circ}}{dt} + \frac{d V_{Go}^{\LARGE \circ}}{dt} = \frac{1}{B_g} q_g + \frac{R_s}{B_o} q_o

volumetric gas-component flow rate in wellbore draining points recalculated to standard surface conditions

q_w = \frac{d V_w}{dt}

volumetric water-phase flow rate in wellbore draining points

q_o = \frac{d V_o}{dt}

volumetric oil-phase flow rate in wellbore draining points

q_g = \frac{d V_g}{dt}

volumetric gas-phase flow rate in wellbore draining points

q^S_W =\frac{dV_{Ww}^S}{dt}

total well volumetric water-component flow rate

q^S_O = \frac{d (V_{Oo}^S + V_{Og}^S )}{dt}

total well volumetric oil-component flow rate

q^S_G = \frac{d (V_{Gg}^S + V_{Go}^S )}{dt}

total well volumetric gas-component flow rate

q^S_L = q^S_W + q^S_O

total well volumetric liquid-component flow rate

P_w = P_w (t, \vec r)

water-phase pressure pressure distribution and dynamics

P_o = P_o (t, \vec r)

oil-phase pressure pressure distribution and dynamics

P_g = P_g (t, \vec r)

gas-phase pressure pressure distribution and dynamics

\vec u_w = \vec u_w (t, \vec r)

water-phase flow speed distribution and dynamics

\vec u_o = \vec u_o (t, \vec r)

oil-phase flow speed distribution and dynamics

\vec u_g = \vec u_g (t, \vec r)

gas-phase flow speed distribution and dynamics

P_{cow} = P_{cow} (s_w)

capillary pressure at the oil-water phase contact as function of water saturation


P_{cog} = P_{cog} (s_ g)

capillary pressure at the oil-gas phase contact as function of gas saturation

k_{rw} = k_{rw}(s_w, \ s_g)

relative formation permeability to water flow as function of water and gas saturation

k_{ro} = k_{ro}(s_w, \ s_g)

relative formation permeability to oil flow as function of water and gas saturation

k_{rg} = k_{rg}(s_w, \ s_g)

relative formation permeability to gas flow as function of water and gas saturation

\phi = \phi(P)

porosity as function of formation pressure

k_a = k_a(P)

absolute formation permeability to air

\vec g = (0, \ 0, \ g)

gravitational acceleration vector

g = 9.81 \ \rm m/s^2

gravitational acceleration constant

\rho_\alpha(P,T)

mass density of \alpha-phase fluid

\mu_\alpha(P,T)

viscosity of \alpha-phase fluid

\lambda_t(P,T,s_w, s_o, s_g)

effective thermal conductivity of the rocks with account for multiphase fluid saturation

\lambda_r(P,T)

rock matrix thermal conductivity

\lambda_\alpha(P,T)

thermal conductivity of \alpha-phase fluid

\rho_r(P,T)

rock matrix mass density

\eta_{s \alpha}(P,T)

differential adiabatic coefficient of \alpha-phase fluid

c_{pr}(P,T)

specific isobaric heat capacity of the rock matrix

c_{p\alpha}(P,T)

specific isobaric heat capacity of \alpha-phase fluid

\epsilon_\alpha (P, T)

differential Joule–Thomson coefficient of \alpha-phase fluid

дифференциальный коэффициент Джоуля-Томсона фазы  \alpha

  • No labels