This paper presents the full set of Volatile Oil equations.

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

Equations

The Volatile Oil flow dynamics is defined by the following set of equations.

\partial_t \bigg [  \phi \ \rho_W  \bigg ]  + \nabla  \bigg ( \rho_{Ww} \ \mathbf{u}_w     \bigg )      = q_{mW}(\mathbf{r})

\partial_t \bigg [  \phi \ \rho_O \bigg ]  + \nabla  \bigg (    \rho_{Oo} \ \mathbf{u}_o 

+   \rho_{Og} \  \mathbf{u}_g \bigg )       = q_{mO}(\mathbf{r})

\partial_t \bigg [  \phi \ \rho_G  \bigg ]  +  \nabla \bigg (   \rho_{Go} \   \mathbf{u}_o

+     \rho_{Gg} \ \mathbf{u}_g  \bigg )     = q_{mG}(\mathbf{r})

\mathbf{u}_w = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ ( \nabla P_w - \rho_w \mathbf{g} )

\mathbf{u}_o = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ (  \nabla P_o - \rho_o \mathbf{g} )

\mathbf{u}_g = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g - \rho_g \mathbf{g} )

P_o - P_w = P_{cow}(s_w)

P_o - P_g = P_{cog}(s_g)

s_w + s_o + s_g = 1

(\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}

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

Equations  –  define the continuity of the fluid components flow or equivalently represent the mass conservation of each mass component  during its transportation in space. 

Equations  –  define the motion dynamics of each phase, represnted as linear correlation between phase flow speed   and partial pressure gradient of this phase  (which is also called Darci flow with account of the gravity and relative permeability).

Equations  –  define the hydrodynamic inter-facial balance between the phases with account of capillary pressure in porous formation . The key assumption is that capillary pressure at oil-water boundary is a function of  water saturation alone  and capillary pressure at oil-gas boundary is a function of  gas saturation alone 

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   implies that porous space is fully occupied by fluid at all times .

Equation   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  defines the speed of change of  heat energy  volumetric density.

In impermeable rocks () heat flow is defined by heat conduction only:

 \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:

(\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:

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

The term  represents heat convection defined by the mass flow. 

The term  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  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  –  represent the system of 16 scalar equations on 16 unknowns: 

,

which are all functions of time and space coordinates .

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

\partial_t \bigg [  \phi \ \rho_W  \bigg ]  + \nabla \bigg (     \rho_w \ \mathbf{u}_w     \bigg )      =  q_{mW}(\mathbf{r}) 

\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})

\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}) 

\mathbf{u}_w = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ ( \nabla P_w - \rho_w  \mathbf{g} )

\mathbf{u}_o = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ ( \nabla P_o - \rho_o   \mathbf{g} )

\mathbf{u}_g = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g - \rho_g  \mathbf{g} )

P_o - P_w = P_{cow}(s_w)

P_o - P_g = P_{cog}(s_g)

s_w + s_o + s_g = 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})

\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})

\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})

\mathbf{u}_w = - k_a \ \frac{k_{rw}(s_w, s_g)}{\mu_w} \ ( \nabla  P_w - \rho_w \  \mathbf{g} )

\mathbf{u}_o = - k_a \ \frac{k_{ro}(s_w, s_g)}{\mu_o} \ ( \nabla P_o - \rho_o \ \mathbf{g} )

\mathbf{u}_g = - k_a \ \frac{k_{rg}(s_w, s_g)}{\mu_g} \ ( \nabla P_g - \rho_g \ \mathbf{g} )

P_o - P_w = P_{cow}(s_w)

P_o - P_g = P_{cog}(s_g)

s_w + s_o + s_g = 1

В уравнениях  –  правые части равны нулю во всем объеме пласта за исключением контакта скважин с пластом, который описывается моделью скважины (см. ниже).

Начальное условие

Начальное условие по температуре задается распределением температурного поля:

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

Начальное условие на давления, скоростей и насыщенности задается одним из двух вариантов.

Условие I  – Стационарный старт

Стационарный старт означает, что до начального момента времени поле давлений  , скоростей   и насыщенностей  находилось в стационарном (не меняющемся во времени) состоянии, соответствующем гидродинамическому равновесию:

\nabla \cdot \bigg (     \frac{1}{B_w} \ \mathbf{u}_w     \bigg )_{t=0}      = 0

\nabla \cdot \bigg (     \frac{1}{B_o} \ \mathbf{u}_o 

+    \frac{R_v}{B_g} \   \mathbf{u}_g    \bigg )_{t=0}      = 0

\nabla \cdot \bigg (     \frac{1}{B_g} \ \mathbf{u}_g

+    \frac{R_s}{B_o} \   \mathbf{u}_o     \bigg )_{t=0}    = 0

\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

\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

\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

P_o(0, \mathbf{r}) - P_w(0, \mathbf{r}) = P_{cow}(s_w)

P_o(0, \mathbf{r}) - P_g(0, \mathbf{r}) = P_{cog}(s_g)

s_w + s_o + s_g = 1

Условие II – Нестационарный старт

Нестационарный старт означает, что к начальному моменту времени поле насыщенностей  является произвольным, с условием

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

поле давлений  является произвольным с условием

P_o(0, \mathbf{r}) - P_w(0, \mathbf{r}) = P_{cow}(s_w)

P_o(0, \mathbf{r}) - P_g(0, \mathbf{r}) = P_{cog}(s_g).

При этом начальное поле скоростей  автоматически рассчитывается по следующим формулам:

\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} )

\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} )

\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} )

На практике, нестационарное начальное поле давлений, скоростей и насыщенностей является, как правило, результатом промежуточных расчетов этой же модели, либо более крупной модели.

Краевое условие на внешней границе

Краевое условие на температурное поле на внешней границе задается одним из двух вариантов

Термодинамическое Условие I  – Фиксированная температура 

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

Термодинамическое Условие II  – Фиксированный теплообмен

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

где   – коэффициент теплообмена на границе резервуара.

Краевое условие на поле давления, скоростей  насыщенности на внешней границе задается одним из двух вариантов

Гидродинамическое Условие I  – Непроницаемая граница 

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

где   – вектор нормали к границе  и .

Гидродинамическое Условие II – Постоянное давление на границе 

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

где  .

Краевое условие на разломах

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

...

Моделирование скважины

Модель притока (или закачки) на каждой скважине связывает объемы добычи (закачки) каждой фазы и перепад давления в пласте и на забое скважины и задается следующей формулой:

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}}

где  – линия контакта скважины и порового коллектора.

Забойное давление  напротив точки контакта скважины и пласта (на глубине )   определяется по одному из трех популярных на практике условий:

• Условие   I – Контроль по забойному давлению
• Условие  II – Контроль по жидкости
• Условие III – Контроль по нефти

Условие I   – Контроль по забойному давлению

Это условие предполагает, что в каждый момент времени известно опорное забойное давление  на глубине  , а забойное давление в каждой точке контакта скважины и пласта  рассчитывается по формуле:

P_{wf}(t) = P_{wf}(t, h_0) +  P_{\delta}(t, \delta h)

где  – изменение забойного давления вдоль ствола скважины в зависимости от характера мультифазного потока в стволе скважины.

При адаптации модели к промысловым данным это условие выполняется для

• нагнетательных скважин, чьи забойные давления пересчитываются по 
  • показаниям устьевых манометров с учетом потерь на трение в стволе скважины
  • по известному давлению на выходе КНС с учетом потерь на трение в стволе скважины и наземных трубопроводах
• для добывающих скважин с мониторингом забойного давления по 
  • глубинным манометрам 
  • эхолотам 
• для добывающих скважин с низким забойным давлением (когда уровень находится вблизи точки подвеса насоса или газлифтного клапана).

При прогнозных расчетах это условие выполняется для

• нагнетательных скважин, чьи режимы закачки определяются давлением на КНС с учетом потерь на трение в стволе скважины и наземных трубопроводах
• для добывающих скважин с низким забойным давлением (когда уровень находится вблизи точки подвеса насоса или газлифтного клапана).

При этом для фонтанной, газлифтной и насосной эксплуатации скважин с забойным давлением выше критического это условие не является физичным и необходимо прогнозировать работу скважины согласно Условию II.

Условие II  – Контроль по жидкости

Это условие предполагает, что известна добыча жидкости на сепараторе каждой скважины  и изменение забойного давления на каждой скважине  рассчитывается по формуле:

P_{wf}(t) = P_{wf}(t, h_{ref}) +  P_{\delta}(t, \delta h)

где  – изменение забойного давления вдоль ствола скважины в зависимости от характера мультифазного потока в стволе скважины,

а опорное давление на глубине   определяется по следующей формуле:

P_{wf}(t, h_{ref}) = 
\frac{ \int_{\Gamma_{WRC}}  \bigg(      
 
\frac{M_w (P_{ew}- \delta P_{wf})}{B^S_w} + \frac{M_o (P_{eo}- \delta P_{wf})}{B^S_o} + \frac{R_v M_g (P_{eg}- \delta P_{wf})}{B^S_g}
 
\bigg)  T_h  dh - q_L(t) }
 
{ \int_{\Gamma_{WRC}}  \bigg(  
 
\frac{M_w}{B^S_w} + \frac{M_o}{B^S_o} + \frac{R_v M_g}{B^S_g}
 
      \bigg)  T_h dh }

которая обеспечивает устьевой дебит по жидкости в размере :

q_L(t) = q_W(t) + q_O(t) = 
 \int_{\Gamma_{WRC}}  \bigg(      

\frac{M_w (P_{ew} - P_{wf}(t, h_{ref}) -  P_{\delta})}{B^S_w} 
+ \frac{M_o (P_{eo} - P_{wf}(t, h_{ref}) - P_{\delta})}{B^S_o} 
+ \frac{R_v M_g (P_{eg} - P_{wf}(t, h_{ref}) - P_{\delta})}{B^S_g}

\bigg)  T_h  dh  

Если пользователь ввел ограничение на минимальное забойное давление  (определяемое например, глубиной спуска насоса или газлифтного клапана), то при достижении  скважина автоматически переходит в режим постоянного давления: 

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}^{min} = const

которое будет сопровождаться изменением дебита всех фаз согласно  – .

Этот режим соотвествует работе насоса с пониженным КПД и в случае если условия на границе контакта поменяются (например, в процессе подъема пластового давления) и потенциал забойного давления согласно  поднимется выше , то скважина опять перейдет в режим работы с заданным дебитом .

При адаптации модели к промысловым данным это условие выполняется для всех типов добывающих и нагнетательных скважин, для которых отборы известны точно (что, кстати, далеко не всегда имеет место быть на практике).

При прогнозных расчетах это условие выполняется для

• нагнетательных и добывающих скважин, чьи режимы работы определяются диаметром штуцером
• для добывающих скважин с высоким забойным давлением (когда уровень находится выше точки подвеса насоса или газлифтного клапана).

При этом для скважин с низким забойным давлением это условие не является физичным и необходимо прогнозировать работу скважины согласно Условию I.

Условие III – Контроль по нефти

Это условие предполагает, что добыча воды и газа неизвестна (или известна неточно) и забойное давление добывающей скважины  в каждый момент времени определяется только значениями устьевых отборов нефти  (которые, как правило, известны точно) по следующей формуле:

P_{wf}(t) = P_{wf}(t, h_{ref}) + P_{\delta}(t, \delta h)

где  – изменение забойного давления вдоль ствола скважины в зависимости от характера мультифазного потока в стволе скважины,

а опорное давление на глубине   определяется по следующей формуле:

P_{wf}(t, h_{ref}) = 
\frac{ \int_{\Gamma_{WRC}}  \bigg(      

 \frac{M_o (P_{eo}- P_{\delta})}{B^S_o} + \frac{R_v M_g (P_{eg}- P_{\delta})}{B^S_g}

\bigg)  T_h  dh - q_O(t) }

{ \int_{\Gamma_{WRC}}  \bigg(  

 \frac{M_o}{B^S_o} + \frac{R_v M_g}{B^S_g}

      \bigg)  T_h dh }

которая обеспечивает устьевой дебит по жидкости в размере :

q_O(t) = 
 \int_{\Gamma_{WRC}}  \bigg(      

\frac{M_o (P_{eo} - P_{wf}(t, h_{ref}) - P_{\delta})}{B^S_o} 
+ \frac{R_v M_g (P_{eg} - P_{wf}(t, h_{ref}) - P_{\delta})}{B^S_g}

\bigg)  T_h  dh 

Если пользователь ввел ограничение на минимальное забойное давление  (определяемое например, глубиной спуска насоса или газлифтного клапана), то при достижении  скважина автоматически переходит в режим постоянного давления

P_o(t, \vec r) \big|_{\Gamma_{WRC}} = P_{wf}^{min} = const

которое будет сопровождаться изменением дебита по нефти согласно .

Этот режим соотвествует работе насоса с пониженным КПД и в случае если условия на границе контакта поменяются (например, в процессе подъема пластового давления) и потенциал забойного давления согласно  поднимется выше , то скважина опять перейдет в режим работы с заданным дебитом .

Условие III по своему определению накладывается только на добывающие скважины и выполняется для всех типов добывающих скважин, для которых отборы нефти известны точно (что наиболее часто встречается на практике).

При этом условие на нагнетательных скважинах не оговорено и может быть как I-ого так и II-ого типа в зависимости от реализации системы ППД.

При прогнозных расчетах условие III использоваться не может в силу своей нефизичности, за исключением случая безводной эксплуатации недонасыщенной нефти, в котором это условие становится физичным и эквивалетным Условию II (контроль по жидкости). 

На практике Условие III рекомендуется накладывать для первичной настройки модели (настройки ее базовых параметров) и потом рекомендуется переключать контроль на Условие I или Условие II в зависимости от промысловых условий эксплуатации скважин.

The list of dynamic flow properties and model parameters