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 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 fo the properties in the above equation is brought in The list of dynamic flow properties and model parameters.
Equations (1) – (3) suggest no sources of flow in the right side except the contacts between wells and reservroir which is specified by well models as boundary conditions (see below).
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 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 and phase velocities \{ \mathbf{u}_w, \ \mathbf{u}_o, \ \mathbf{u}_g \} were zero, corresponding to hydrodynamic equilibrium:
(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
Equilibrium Start means that phase pressure \{ P_w, \ P_o, \ P_g \}, phase velocities \{ \mathbf{u}_w, \ \mathbf{u}_o, \ \mathbf{u}_g \} and phase saturations \{ s_w, \ s_o, \ s_g, \} were in stationary (not varying in time) conditions, corresponding to hydrodynamic equilibrium:
Нестационарный старт означает, что к начальному моменту времени поле насыщенностей \{ s_w, \ s_o, \ s_g, \} является произвольным, с условием
(42) | s_w(0, \mathbf{r}) + s_o(0, \mathbf{r}) + s_g(0, \mathbf{r}) = 1 |
поле давлений \{ P_w, \ P_o, \ P_g \} является произвольным с условием
(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). |
При этом начальное поле скоростей \{ \vec u_w, \ \vec u_o, \ \vec u_g \} автоматически рассчитывается по следующим формулам:
(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} ) |
На практике, нестационарное начальное поле давлений, скоростей и насыщенностей является, как правило, результатом промежуточных расчетов этой же модели, либо более крупной модели.
Краевое условие на внешней границе
Краевое условие на температурное поле на внешней границе задается одним из двух вариантов
Термодинамическое Условие I – Фиксированная температура
T(t, \mathbf{r}) |_{\Gamma_e} = T_e( \mathbf{r}) |
Термодинамическое Условие II – Фиксированный теплообмен
(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) |
где \zeta – коэффициент теплообмена на границе резервуара.
Краевое условие на поле давления, скоростей насыщенности на внешней границе задается одним из двух вариантов
Гидродинамическое Условие I – Непроницаемая граница \Gamma_e
(49) | \big( \mathbf{n}, \ (\nabla P_\alpha(t, \mathbf{r}) - \rho_\alpha \mathbf{r}) \big) \big|_{\Gamma_e} = 0 |
где \mathbf{n} – вектор нормали к границе \Gamma_e и \alpha = \{ w, o, g \}.
Гидродинамическое Условие II – Постоянное давление на границе \Gamma_e
(50) | P_\alpha(t, \mathbf{r}) \big|_{\Gamma_e} = P_i = const |
где \alpha = \{ w, o, g \}.
Краевое условие на разломах
Предполагается выполнение одного из двух условий на разломах (индивидуально по каждому разлому).
Условие I – Непроницаемый разлом \Gamma_F
(51) | \big( \mathbf{n}, \ ( \nabla P_\alpha(t, \mathbf{r}) - \rho_\alpha \mathbf{g}) \big) \big|_{\Gamma_F} = 0 |
где \vec n – вектор нормали к разлому \Gamma_F и \alpha = \{ w, o, g \}.
Условие II – Проницаемый разлом \Gamma_F
(52) | ... |
где \alpha = \{ w, o, g \}.
Моделирование скважины
Модель притока (или закачки) на каждой скважине связывает объемы добычи (закачки) каждой фазы и перепад давления в пласте и на забое скважины и задается следующей формулой:
(53) | 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}} |
где \Gamma_{WRC} – линия контакта скважины и порового коллектора.
Забойное давление P_{wf}(t) \big|_{\Gamma_{WFC}} = P_{wf}(t,h) напротив точки контакта скважины и пласта (на глубине h) определяется по одному из трех популярных на практике условий:
- Условие I – Контроль по забойному давлению
- Условие II – Контроль по жидкости
- Условие III – Контроль по нефти
Условие I – Контроль по забойному давлению
Это условие предполагает, что в каждый момент времени известно опорное забойное давление P_{wf} (t, h_{ref}) на глубине h_{ref} , а забойное давление в каждой точке контакта скважины и пласта рассчитывается по формуле:
(54) | P_{wf}(t) = P_{wf}(t, h_0) + P_{\delta}(t, \delta h) |
где P_{\delta}(t, \delta h) – изменение забойного давления вдоль ствола скважины в зависимости от характера мультифазного потока в стволе скважины.
При адаптации модели к промысловым данным это условие выполняется для
- нагнетательных скважин, чьи забойные давления пересчитываются по
- показаниям устьевых манометров с учетом потерь на трение в стволе скважины
- по известному давлению на выходе КНС с учетом потерь на трение в стволе скважины и наземных трубопроводах
- для добывающих скважин с мониторингом забойного давления по
- глубинным манометрам
- эхолотам
- для добывающих скважин с низким забойным давлением (когда уровень находится вблизи точки подвеса насоса или газлифтного клапана).
При прогнозных расчетах это условие выполняется для
- нагнетательных скважин, чьи режимы закачки определяются давлением на КНС с учетом потерь на трение в стволе скважины и наземных трубопроводах
- для добывающих скважин с низким забойным давлением (когда уровень находится вблизи точки подвеса насоса или газлифтного клапана).
При этом для фонтанной, газлифтной и насосной эксплуатации скважин с забойным давлением выше критического это условие не является физичным и необходимо прогнозировать работу скважины согласно Условию II.
Условие II – Контроль по жидкости
Это условие предполагает, что известна добыча жидкости на сепараторе каждой скважины q_L(t) и изменение забойного давления на каждой скважине P_{wf} (t) рассчитывается по формуле:
(55) | P_{wf}(t) = P_{wf}(t, h_{ref}) + P_{\delta}(t, \delta h) |
где P_{\delta}(t, \delta h) – изменение забойного давления вдоль ствола скважины в зависимости от характера мультифазного потока в стволе скважины,
а опорное давление на глубине h_{ref} определяется по следующей формуле:
(56) | 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):
(57) | 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_{wf}^{ \ min} (определяемое например, глубиной спуска насоса или газлифтного клапана), то при достижении P_{wf}(t) = P_{wf}^{ \ min} скважина автоматически переходит в режим постоянного давления:
(58) | 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 |
которое будет сопровождаться изменением дебита всех фаз согласно
– .Этот режим соотвествует работе насоса с пониженным КПД и в случае если условия на границе контакта поменяются (например, в процессе подъема пластового давления) и потенциал забойного давления согласно (55) поднимется выше P_{wf}^{min}, то скважина опять перейдет в режим работы с заданным дебитом q_L(t).
При адаптации модели к промысловым данным это условие выполняется для всех типов добывающих и нагнетательных скважин, для которых отборы известны точно (что, кстати, далеко не всегда имеет место быть на практике).
При прогнозных расчетах это условие выполняется для
- нагнетательных и добывающих скважин, чьи режимы работы определяются диаметром штуцером
- для добывающих скважин с высоким забойным давлением (когда уровень находится выше точки подвеса насоса или газлифтного клапана).
При этом для скважин с низким забойным давлением это условие не является физичным и необходимо прогнозировать работу скважины согласно Условию I.
Условие III – Контроль по нефти
Это условие предполагает, что добыча воды и газа неизвестна (или известна неточно) и забойное давление добывающей скважины P_{wf}(t,h) в каждый момент времени определяется только значениями устьевых отборов нефти q_O(t) (которые, как правило, известны точно) по следующей формуле:
(59) | P_{wf}(t) = P_{wf}(t, h_{ref}) + P_{\delta}(t, \delta h) |
где – изменение забойного давления вдоль ствола скважины в зависимости от характера мультифазного потока в стволе скважины,
а опорное давление на глубине h_{ref} определяется по следующей формуле:
(60) | 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):
(61) | 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_{wf}^{ \ min} (определяемое например, глубиной спуска насоса или газлифтного клапана), то при достижении P_{wf}(t) = P_{wf}^{ \ min} скважина автоматически переходит в режим постоянного давления
(62) | P_o(t, \vec r) \big|_{\Gamma_{WRC}} = P_{wf}^{min} = const |
которое будет сопровождаться изменением дебита по нефти согласно (61).
Этот режим соотвествует работе насоса с пониженным КПД и в случае если условия на границе контакта поменяются (например, в процессе подъема пластового давления) и потенциал забойного давления согласно (59) поднимется выше P_{wf}^{min}, то скважина опять перейдет в режим работы с заданным дебитом q_O(t).
Условие III по своему определению накладывается только на добывающие скважины и выполняется для всех типов добывающих скважин, для которых отборы нефти известны точно (что наиболее часто встречается на практике).
При этом условие на нагнетательных скважинах не оговорено и может быть как I-ого так и II-ого типа в зависимости от реализации системы ППД.
При прогнозных расчетах условие III использоваться не может в силу своей нефизичности, за исключением случая безводной эксплуатации недонасыщенной нефти, в котором это условие становится физичным и эквивалетным Условию II (контроль по жидкости).
На практике Условие III рекомендуется накладывать для первичной настройки модели (настройки ее базовых параметров) и потом рекомендуется переключать контроль на Условие I или Условие II в зависимости от промысловых условий эксплуатации скважин.
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 \ 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 |