Displacement of hydrocarbon liquids by water using the models of zonally heterogeneous deformable formations

Толық мәтін

Введение

Известно, что разработка глубокозалегающих газоконденсатных и нефтяных залежей сопровождается деформацией пород коллекторов, в результате чего изменяются их емкостные и коллекторские характеристики [1, 2]. Установлено, что при более широком диапазоне изменения пластового давления деформация горных пород может иметь существенно нелинейный характер [3–5]. Кроме этого, при этом могут проявляться ползучесть горных пород [6–8]. Причем в одном и том же пласте в зависимости от значения внутрипорового давления деформации скелета коллектора могут показать себя по-разному [9–11]. Так, если вблизи призабойной зоны, где пластовое давление намного ниже его первоначального значения, деформации пласта-коллектора происходят по одному закону, а на контуре (или в дали от забоя), где сравнительно высокое давление (или давление выше определенного предела) скелет коллектора сжимается по другому закону. Помимо этого, иногда пласт имеет по проницаемости зональную неоднородность. В работах [12, 13] получены решения задач моделирования процессов истощения летучих нефтей и газоконденсатных залежей на зонально-неоднородных коллекторах. При этом пласт состоял из двух зон, отличающихся коллекторскими свойствами. Данная работа является продолжением исследований, начатых в отмеченных работах. Она посвящена более сложной задаче – исследованию процесса вытеснения к эксплуатационной скважине углеводородных систем в зонально-неоднородных пластах. При этом пласт круговой формы, разрабатываемый одной центральной скважиной, состоит из двух зон с различными коллекторско-емкостными и реологическими свойствами. Учитывается неполнота вытеснения, сжимаемость воды и PVT свойства углеводородной системы – фазовое превращение, массообмен между фазами углеводородной системы, что отличает предложенное решение от существующих.

Объекты и методика исследования

Схематическое изображение течения воды при вытеснении углеводородов в зонально-неоднородном пласте иллюстрируется на рис. 1. Согласно изображению, обводненная часть пласта состоит из двух зон – зоны c проницаемостью k₁ (зона I) и зоны c проницаемостью k₂ (зона II), которые первоначально были насыщенными нефтью. К рассматриваемому моменту текущее положение водонефтяного контакта (ВНК) r_v [14–16].

Уравнение нестационарной фильтрации воды в зоне I, т. е. в области r_k£r£Rr_k, будет иметь следующий вид [17] (1):

$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial P_1}{\partial r}\right)=\frac{1}{a}\frac{\partial P_1}{\partial t},$  (1)

где

$\begin{array}{c}{P}_1=\int k_1{\text{γ}}_{1}d{p}_1,a=\frac{k_1\left(P_0\right)}{\left({\text{β}}_k+{\text{β}}_v\right){\mu }_v{m}_1\left(P_0\right)B^{*}},\\ P_2=\int k_2{\gamma }_{1}d{p}_2,{\beta }_k, {\beta }_v\end{array}$

– соответственно, коэффициенты изменения проницаемости и сжимаемости воды, m_v – вязкость воды, m₁ – пористость пласта в зоне I, g₁ – плотность воды. B^* является угловым коэффициентом ломанного участка кривой зависимости $\overline{m}\overline{\gamma }=\overline{m\gamma }\left(P_1\right).$

Уравнение (1) решается при следующих начальных и граничных условиях (2)–(4):

$P_1\left(r\mathrm{,0}\right)=P_0,$ (2)

$P_1\left(R_k,t\right)=P_k\left(t\right),$ (3)

$P_1\left(r_k,t\right)=P^{\text{'}}\left(t\right).$ (4)

Течение воды в зоне I, согласно результатам предыдущих исследований, можем принимать стационарным, что позволит упростить задачу. При этом уравнение фильтрации воды в зоне I, т. е. в области r_v£r£r_k, напишем в виде (5)

$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial P_2}{\partial r}\right)=\mathrm{0,}$ (5)

при следующих граничных условиях (6), (7):

$P_2\left(r_k,t\right)=P^{\text{'}}\left(t\right),$ (6)

$P_2\left(r_v,t\right)=P_v\left(t\right).$ (7)

 

Рис. 1. Схематическое изображение течения воды при вытеснении углеводородной системы в зонально-неоднородном пласте

Fig. 1. Schematic representation of water flow when displacing a hydrocarbon system in a zonally heterogeneous formation

 

Имеем дополнительное условие на R_k (8)

${\left.\frac{\partial P_1}{\partial r}\right|}_{r=R_k}=\frac{q_{vz}}{A\left(P_k\right)},$

где

$A\left(P_k\right)=\frac{2\pi R_{k}h{F}_v{k}_0}{\left({\beta }_k+{\beta }_v\right){\mu }_v{e}^{{\beta }_v(p_k-p_{0)}}},$ (8)

и условие неразрывности (9):

${\left.\frac{\partial P_1}{\partial r}\right|}_{r=r_k}={\left.\frac{\partial P_2}{\partial r}\right|}_{r=r_k},$ (9)

где R_k – радиус контура заводнения; h – толщина контура; F_v – относительная площадь текущего положения водонефтяного контакта; k₀ – проницаемость на нулевом контуре.

Решением системы уравнений (1) и (5) при краевых условиях (2)–(4), (6)–(9) получено следующее выражение для определения мгновенного значения расхода вторгшейся в продуктивную часть пласта воды (10):

$q_v=\frac{2\pi R_{k}h{F}_v{k}_0}{\left({\beta }_k+{\beta }_v\right){\mu }_v{e}^{{\beta }_v(p_k-p_{0)}}}\frac{P^{\text{'}}-P_v}{\ln \frac{r_k}{r_v}}.$ (10)

Здесь P_v будет определяться ниже решением задачи фильтрации нефти (или газа) к скважине по пластовому давлению, а P' – значение функции P₁ на r_k – определяется по следующему выражению, полученному с помощью условия (9), (11):

$\begin{array}{c}{\overline{P}}^{\text{'}}=\frac{\left(\begin{array}{l}-2{\displaystyle \sum }_{n=1}^{\infty }\frac{J_0\left({\overline{R}}_k{\alpha }_n\right)e^{-{a^2}_{n}t}}{J_0^2\left({\overline{R}}_{n0}{\alpha }_n\right)-J_0^2\left({\overline{R}}_k{\alpha }_n\right)}\times \\ \times \left\{\begin{array}{l}\{{\overline{P}}_0-P_k\}\left[J_0\right({\overline{R}}_{n0}{\alpha }_n)-J_0({\overline{R}}_k{\alpha }_n\left)\right]-\\ -{\overline{P}}_0{J}_0\left({\overline{R}}_k{\alpha }_n\right)\end{array}\right\}\end{array}\right)}{2{\displaystyle \sum _n^{\infty }\frac{J_0^2\left({\overline{R}}_k{\alpha }_n\right)e^{-a_n^{2}t}}{J_0^2\left({\overline{R}}_{n0}{\alpha }_n\right)-J_0^2\left({\overline{R}}_k{\alpha }_n\right)}+\frac{1}{\ln \frac{{\overline{R}}_k}{{\overline{R}}_{n0}}}+\frac{F_v}{\ln \frac{{\overline{R}}_{n0}}{{\overline{R}}_n}}}}\to \\ \to \frac{\frac{{\overline{P}}_k}{\ln \frac{{\overline{R}}_k}{{\overline{R}}_{n0}}}+\frac{{\overline{P}}_n}{\ln \frac{{\overline{R}}_{n0}}{{\overline{R}}_n}}}{2{\displaystyle \sum _{n=1}^{\infty }\frac{J_0^2\left({\overline{R}}_k{\alpha }_n\right)e^{-a_n^{2}t}}{J_0^2\left({\overline{R}}_{n0}{\alpha }_n\right)-J_0^2\left({\overline{R}}_k{\alpha }_n\right)}+\frac{1}{\ln \frac{{\overline{R}}_k}{{\overline{R}}_{n0}}}+\frac{F_v}{\ln \frac{{\overline{R}}_{n0}}{{\overline{R}}_n}}}},\end{array}$

где J₀ – момент инерции рассматриваемого сектора зоны; a_n – угол отклонения контура от вертикали; ${\stackrel{\_\_}{P}}_0$ – давление в начале контура; ${\stackrel{\_\_}{P}}_k$ – давление в конце контура.

Выражение для определения значения P₁ на контуре закачки P_k получено в виде:

$\begin{array}{c}{\overline{P}}_k=\frac{\left[\begin{array}{l}\frac{q_{vz}}{A\left(P_k\right)}+\frac{{\overline{P}}^{\text{'}}}{\ln \frac{{\overline{R}}_k}{{\overline{R}}_{n0}}}\frac{1}{{\overline{R}}_k}+\\ +\frac{2}{{\overline{R}}_k}{\displaystyle \sum }_{n=1}^{\infty }\frac{J_0\left({\overline{R}}_k{\alpha }_n\right)e^{-a_n^{2}t}}{J_0^2\left({\overline{R}}_{n0}{\alpha }_n\right)-J_0^2\left({\overline{R}}_k{\alpha }_n\right)}\end{array}\right]}{\frac{1}{\ln \frac{{\overline{R}}_k}{{\overline{R}}_{n0}}}\frac{1}{{\overline{R}}_k}+\frac{2}{{\overline{R}}_k}{\displaystyle \sum }_{n=1}^{\infty }\frac{J_0\left({\overline{R}}_k{\alpha }_n\right)e^{-a_n^{2}t}}{J_0^{}\left({\overline{R}}_{n0}{\alpha }_n\right)-J_0^{}\left({\overline{R}}_k{\alpha }_n\right)}}\to \\ \to \frac{\left\{\begin{array}{l}{\overline{P}}_0\left[J_0\right({\overline{R}}_{n0}{\alpha }_n)-J_0({\overline{R}}_k{\alpha }_n\left)\right]-\\ -({\overline{P}}_0-{\overline{P}}^{\text{'}})J_0\left({\overline{R}}_k{\alpha }_n\right)\end{array}\right\}}{\frac{1}{\ln \frac{{\overline{R}}_k}{{\overline{R}}_{n0}}}\frac{1}{{\overline{R}}_k}+\frac{2}{{\overline{R}}_k}{\displaystyle \sum }_{n=1}^{\infty }\frac{J_0\left({\overline{R}}_k{\alpha }_n\right)e^{-a_n^{2}t}}{J_0^{}\left({\overline{R}}_{n0}{\alpha }_n\right)-J_0^{}\left({\overline{R}}_k{\alpha }_n\right)}},\end{array}$

где переход от функции P к давлению p осуществляется следующим выражением

$p=p_0+\frac{\ln P}{{\beta }_k+{\beta }_v}.$

Для неполноты вытеснения между объемом пор заводненной части $V$ и объемом вторгшейся в продуктивную часть пласта воды V_v существует соотношение $V=\frac{V_v}{1-{\text{ρ}}_{\text{ost}}}.$ Зная мгновенное значение расхода вторгшейся в продуктивную часть пласта воды, вычисленное по (10) с учетом (11), для объема пор заводненной части пласта V и следовательно для текущего радиуса ВНК r_v напишем:

$\begin{array}{c}V=\frac{{\displaystyle {\int }_0^T{q}_{v}dt}}{1-{\text{ρ}}_{\text{ost}}}\approx \frac{{\displaystyle {\sum }_{i=1}^n{q}_{vi}\left(t\right)\Delta t}}{1-{\text{ρ}}_{\text{ost}}},\\ r_v^{2}hm\left(p_v,t\right)\approx \frac{{\displaystyle {\sum }_{i=1}^n{q}_{vi}\left(t\right)\Delta t}}{1-{\text{ρ}}_{\text{ost}}}\\ \text{и\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{r}_v\approx \sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{q}_{vi}\left(t\right)\Delta t}}{\pi hm\left(p_v,t\right)\left(1-{\text{ρ}}_{\text{ost}}\right)}}.\end{array}$

Теперь рассмотрим процесс фильтрации двухфазной углеводородной системы в деформируемом пористом пласте [18–20].

Фильтрация углеводородной системы к центральной скважине при вытеснении водой в зонально-неоднородном пласте схематически показана на рис. 2. По схеме видно, что текущее положение фронта воды имеет радиуса r_v. Контур заводнения имеет радиус R_k, а r_v – текущее положение ВНК [21].

Известно, что движение двухфазных углеводородных систем в деформируемых коллекторах представляется сложными нелинейными дифференциальными уравнениями в частных производных [22–24]. Аналитическое решение таких уравнений возможно лишь с применением особых подходов [25]. В настоящей работе для линеаризации уравнений будем применять метод осреднения и функцию Христиановича аналогично работе [26] и ниже получим решение задачи фильтрации двухфазных углеводородных систем к скважине в зонально-неоднородном пористом пласте при вытеснении углеводородной системы водой.

И так, радиальное движение углеводородной системы в зоне II, т. е. в области r_s£r£r_k, описываются уравнениями (12)–(14):

$\frac{1}{r}\frac{\partial }{\partial r}\left[r\varphi \left(p,s\right)\frac{\partial p}{\partial r}\right]=-\frac{\partial }{\partial t}f\left(p,s\right),$ (12)

$\frac{1}{r}\frac{\partial }{\partial r}\left[r{\varphi }_g\left(p,s\right)\frac{\partial p}{\partial r}\right]=-\frac{\partial }{\partial t}{f}_g\left(p,s\right),$ (13)

где

$\begin{array}{c}\varphi \left(p,s\right)=\left[\frac{k_{ro}\left(s\right)}{{\mu }_o\left(p\right)B_o\left(p\right)}+\frac{k_{rg}\left(s\right)p\beta c\left(p\right)}{{\mu }_g\left(p\right)z\left(p\right)p_{at}}\right]k\left(p\right),\\ {\varphi }_g\left(p,s\right)=\left[\begin{array}{l}\frac{k_{rg}\left(s\right)p\beta \left[1-c\left(p\right)\overline{\gamma }\left(p\right)\right]}{{\mu }_g\left(p\right)z\left(p\right)p_{at}}+\\ +\frac{k_{ro}\left(s\right)S\left(p\right)}{{\mu }_o\left(p\right)B_o\left(p\right)}\end{array}\right]k\left(p\right),\\ f\left(p,s\right)=\left[\frac{s}{B_o\left(p\right)}+\left(1-s\right)\frac{p\beta c\left(p\right)}{z\left(p\right)p_{at}}\right]\phi \left(p\right),\\ f_g\left(p,s\right)=\left[\frac{\left(1-s\right)p\beta \left[1-c\left(p\right)\overline{\gamma }\left(p\right)\right]}{z\left(p\right)p_{at}}+s\frac{S\left(p\right)}{B_o\left(p\right)}\right]\phi \left(p\right),\end{array}$  (14)

где k_ro(s), k_rп(s) – относительные фазовые проницаемости для жидкой фазы (например, для нефти или жидкого конденсата в случае фильтрации газоконденсатной смеси) и газовой фазы, соответственно; s – насыщенность пор жидкой фазой (нефтью или конденсатом); z, b – коэффициенты сверхсжимаемости и температурной поправки для газовой фазы; с – содержание потенциально жидких углеводородов в газовой фазе; m_o, m_g – вязкости жидкой и газовой фаз, соответственно; B_o – объемный коэффициент жидкой фазы; S – количество растворенного газа в жидкой фазе; $\overline{\gamma }=\frac{{\gamma }_o\left(p\right)}{{\gamma }_g\left(p\right)}$ – отношение удельных весов жидкой и газовой фаз при пластовом давлении ; p_at – атмосферное давление; j и j₁ – текущее значение эффективной пористости II и I зон пласта, соответственно; k и k₁ – текущее значение эффективной проницаемости II и I зон пласта, соответственно; r – радиальная координата; t – время.

Уравнение (12) описывает нестационарное движение жидких углеводородов и потенциального конденсата в газовой фазе, а (13) – движение газа и паров более легких углеводородов в пористой среде [27]. Аналогичные уравнения движения углеводородной системы в зоне I в области r_k£r£r_v выписываются в следующем виде (15), (16), (16^*):

$\frac{1}{r}\frac{\partial }{\partial r}\left[r\varphi \left(p_1,s_1\right)\frac{\partial p_1}{\partial r}\right]=-\frac{\partial }{\partial t}\left[f\left(p_1,s_1\right)\right],$ (15)

$\frac{1}{r}\frac{\partial }{\partial r}\left[r{\varphi }_g\left(p_1,s_1\right)\frac{\partial p_1}{\partial r}\right]=-\frac{\partial }{\partial t}\left[f_g\left(p_1,s_1\right)\right],$  (16)

$\begin{array}{c}\left(p_1,s_1\right)=\left[\frac{k_{ro}\left(s_1\right)}{{\mu }_o\left(p_1\right)B_o\left(p_1\right)}+\frac{k_{rg}\left(s_1\right)p\beta c\left(p_1\right)}{{\mu }_g\left(p_1\right)z\left(p_1\right)p_{at}}\right]k_1\left(p_1\right),\\ {\varphi }_g\left(p_1,s_1\right)=\left[\begin{array}{l}\frac{k_{rg}\left(s_1\right)p_1\beta \left[1-c\left(p_1\right)\overline{\gamma }\left(p_1\right)\right]}{{\mu }_g\left(p_1\right)z\left(p_1\right)p_{at}}+\\ +\frac{k_{ro}\left(s_1\right)S\left(p_1\right)}{{\mu }_o\left(p_1\right)B_o\left(p_1\right)}\end{array}\right]k_1\left(p_1\right),\end{array}$

$\begin{array}{c}f\left(p_1,s_1\right)=\left[\frac{s_1}{B_o\left(p_1\right)}+\left(1-s_1\right)\frac{p\beta c\left(p_1\right)}{z\left(p_1\right)p_{at}}\right]{\phi }_1\left(p_1\right),\\ f_g\left(p_1,s_1\right)=\left[\begin{array}{l}\frac{\left(1-s_1\right)p_1\beta \left[1-c\left(p_1\right)\overline{\gamma }\left(p_1\right)\right]}{z\left(p_1\right)p_{at}}+\\ +s_1\frac{S\left(p_1\right)}{B_o\left(p_1\right)}\end{array}\right]{\phi }_1\left(p_1\right),\end{array}$ (16*)

где p₁, s₁ – средние давление и насыщенность пор жидкой фазой в зоне I, соответственно; k₁(p₁), j₁(p₁)– эффективные проницаемость и пористость зоны I при давлении p₁.

 

Рис. 2. Схематическое изображение течения углеводородной системы при вытеснении водой в зонально-неоднородном пласте

Fig. 2. Schematic representation of a hydrocarbon system flow when displaced by water in a zonally heterogeneous formation

 

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

Уравнения (12)–(16) являются нелинейными уравнениями, для линеаризации которых, как отмечено выше, применим метод осреднения [28]. Если усреднить пластовое давление по координате r, правые стороны уравнений будут зависеть только от времени [29]. Учитывая это, правую часть уравнений приравниваем некоторой функции Ф(t). Введя функцию, аналогичную функции Хрестиановича уравнения движения в зонах с проницаемостью, соответственно, k₂ и k₁ перепишем в следующем виде (17), (18):

$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial H}{\partial r}\right)=-\text{Ф}\left(t\right),$ (17)

$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial H_1}{\partial r}\right)=-{\text{Ф}}_1\left(t\right),$ (18)

где H, H₁ являются функциями Христиановича H=òf(p,s)dp+const, H₁=òf(p₁,s₁)dp₁+const; Ф(t), Ф₁(t) – некоторые функции, зависящие только от времени и определяемые для фиксированного времени с помощью дополнительных условий.

Системы уравнений (17) и (18) решаются при следующих краевых условиях (19):

$\begin{array}{c}r=r_s,H=H_s;\\ r=r_k,H=H_k;\\ r=r_v,H_1=H_v.\end{array}$ (19)

Дополнительно имеем следующие условия и обозначения (20):

$\begin{array}{c}r=r_k,\frac{\partial H}{\partial r}=\frac{\partial H_1}{\partial r};\\ r=r_v,\frac{\partial H_1}{\partial r}=\frac{q_v}{A_1\left(p_v\right)}.\end{array}$ (20)

Общее решение уравнения (17) при граничных условиях (19) легко получить в виде (21):

$\begin{array}{c}H=\frac{1}{4}\text{Ф}\left(t\right)\left(-r^2+\frac{r_k^2-r_s^2}{\ln \frac{r_k}{r_s}}\ln \frac{r}{r_k}+r_k^2\right)+\\ +\frac{H_k-H_s}{\ln \frac{r_k}{r_s}}\ln \frac{r}{r_k}+H_k.\end{array}$ (21)

Аналогичным образом решается (18), и получается следующее общее решение (22):

$\begin{array}{c}{H}_1=-\frac{1}{4}{\text{Ф}}_1\left(t\right)\left(r^2-r_v^2-\frac{r_v^2-r_k^2}{\ln \frac{r_v}{r_k}}\ln \frac{r}{r_v}\right)+\\ +\frac{H_v-H_k}{\ln \frac{r_v}{r_k}}\ln \frac{r}{r_v}+H_k.\end{array}$  (22)

Из (21) и (22) с учетом условий (20) можно определять Ф(t), Ф₁(t) (23):

$\begin{array}{c}\text{Ф}\left(t\right)=\frac{\left[\begin{array}{l}\frac{\frac{2\left(H_v-H_k\right)}{\ln \frac{r_v}{r_k}}-\frac{q_v{\mu }_v}{\pi h{F}_v}}{r_v^2-\frac{r_v^2-r_k^2}{2\ln \frac{r_v}{r_k}}}\left(r_k^2-\frac{r_v^2-r_k^2}{2\ln \frac{r_v}{r_k}}\right)-\\ -\frac{2\left(H_v-H_k\right)}{\ln \frac{r_v}{r_k}}+\frac{2\left(H_k-H_s\right)}{\ln \frac{r_k}{r_s}}\end{array}\right]}{r_k^2-\frac{r_k^2-r_s^2}{2\ln \frac{r_k}{r_s}}},\\ {\text{Ф}}_1\left(t\right)=\frac{\frac{2\left(H_v-H_k\right)}{\ln \frac{r_v}{r_k}}-\frac{q_v{\mu }_v}{\pi h{F}_v}}{r_v^2-\frac{r_v^2-r_k^2}{2\ln \frac{r_v}{r_k}}}.\end{array}$  (23)

Если учитывать, что дебит скважины $q=2\pi r_{s}h{\left.\frac{\partial H}{\partial r}\right|}_{r-r_s},$ выражение для определения мгновенного дебита скважины получим из (21) с учетом (19) в следующем виде (24):

$q=2\pi h\left[\begin{array}{l}\frac{\left[\begin{array}{l}\frac{H_v-H_k}{r_v^2\ln \frac{r_v}{r_k}-\frac{1}{2}\left(r_v^2-r_k^2\right)}\left(r_k^2-\frac{r_v^2-r_k^2}{2\ln \frac{r_v}{r_k}}\right)-\\ -\frac{H_v-H_k}{\ln \frac{r_v}{r_k}}+\frac{H_k-H_s}{\ln \frac{r_k}{r_s}}\end{array}\right]}{\frac{r_k^2-r_s^2}{2\ln \frac{r_k}{r_s}}-r_k^2}\times \\ \times \left(r_s^2-\frac{r_k^2-r_s^2}{2\ln \frac{r_k}{r_s}}\right)+\frac{H_k-H_s}{\ln \frac{r_k}{r_s}}\end{array}\right].$  (24)

Для перетока из зоны I в зону II на границе r_k получено следующее выражение (25):

$q_1=\pi h\left[2\frac{H_k-H_s}{\ln \frac{r_k}{r_s}}-\text{Ф}\left(t\right)\left(r_k^2-\frac{r_k^2-r_s^2}{2\ln \frac{r_k}{r_s}}\right)\right],$   (25)

где Ф(t) определяется по (23). 

(24) действителен при r_v>r_k. В случае, когда r_v становится равно или меньше r_k, т. е. когда ВНК входит в зону с проницаемостью k₂, дебит скважины определяется в области r_s£r£r_k решением уравнения (17) при граничных условиях (рис. 3) (26):

$\begin{array}{c}r=r_s,H=H_s;\\ r=r_v,H=H_v;\\ r=r_v,\frac{\partial H}{\partial r}=\frac{q_v}{A\left(p_v\right)}.\end{array}$ (26)

 

Рис. 3. Схематическое изображение процесса вытеснения, когда r_v<r_k в области r_s≤r≤r_k

Fig. 3. Schematic representation of displacement when r_v<r_k in the region r_s^≤r≤r_k

 

Выражение для вычисления текущего дебита скважины получено аналогично предыдущему случаю в следующем виде [30] (27):

$q=\pi h\left[\frac{\frac{q_v}{A\left(p_v\right)}-\frac{H_v-H_s}{r_v\ln \frac{r_v}{r_s}}}{\frac{r_v}{2}-\frac{r_v^2-r_s^2}{4r_v\ln \frac{r_v}{r_s}}}\left(r_s^2-\frac{r_v^2-r_s^2}{2\ln \frac{r_v}{r_s}}\right)+2\frac{H_v-H_s}{\ln \frac{r_v}{r_s}}\right].$ (27)

При применении выражений (24), (25) и (27) необходимы определения разницы псевдодавлений H_v–H_k, H_k–H_s и H_v–H_s. Для этой цели применяем аппроксимацию подынтегральной функции  логарифмической функцией вида (28)

$\varphi =a\ln \left(p\right)-b,$ (28)

где коэффициенты a и b находятся из граничных значений функции f по нижеприведенным выражениям. Точность этой аппроксимации подробно исследована в работах [31–33], поэтому не будем уделять внимание.

С учетом этой аппроксимации проинтегрируем функции Христиановича H=òf(p,s)dp+const, H₁=òf(p₁,s₁)dp₁+const в пределах давлений [p_k,p_v], [p_s,p_k], [p_s,p_v] и получим соответствующие выражения для H_v–H_k, H_k–H_s и H_v–H_s в следующем виде (29):

$\begin{array}{c}{H}_v-H_k=a\left(p_v\ln \frac{p_v^{p_v}}{p_k^{p_k}}-p_v+p_k\right)-b\left(p_v-p_k\right),\\ H_k-H_s=a\left(p_k\ln \frac{p_k^{p_k}}{p_s^{p_s}}-p_k+p_s\right)-b\left(p_k-p_s\right),\\ H_v-H_s=a\left(p_v\ln \frac{p_v^{p_v}}{p_s^{p_s}}-p_v+p_s\right)-b\left(p_v-p_s\right),\end{array}$  (29)

где соотношения для вычисления коэффициентов a и b получены из (28) с учетом соответствующих граничных значений f в следующем виде (30):

$\begin{array}{c}a=\frac{{\varphi }_v-{\varphi }_k}{\ln \frac{p_v}{p_k}},b=\frac{{\varphi }_v-{\varphi }_k}{\ln \frac{p_v}{p_k}}\ln p_v-{\varphi }_v,\\ a=\frac{{\varphi }_k-{\varphi }_s}{\ln \frac{p_k}{p_s}},b=\frac{{\varphi }_k-{\varphi }_s}{\ln \frac{p_k}{p_s}}\ln p_k-{\varphi }_s;\\ a=\frac{{\varphi }_v-{\varphi }_s}{\ln \frac{p_v}{p_s}},b=\frac{{\varphi }_v-{\varphi }_s}{\ln \frac{p_v}{p_s}}\ln p_v-{\varphi }_v.\end{array}$   (30)

Здесь f_v, f_k и f_s являются значениями f при давлениях p_v, p_k и p_s, соответственно.

Однако для реализации изложенных соотношений потребуются определения пластовых давлений и насыщенностей пор жидкой фазой на границе между рассматриваемыми зонами и на забое в каждый момент времени. Для этого будем использовать уравнение материального баланса жидкости и газа [34–36].

Определим средние значения временных параметров пластовых давлений и насыщенностей пор жидкой фазой. Полученные выше решения позволяют определять мгновенное значение дебита скважины, т. е. значение, которое соответствует моменту некоторого значения пластовых параметров, таких как пластовое давление, насыщенность пор жидкой фазой, положение ВНК [37]. Для прогнозирования дебита скважины необходим алгоритм для определения отмеченных пластовых параметров в любой момент времени [38–40]. Для этой цели будем использовать уравнения материального баланса для газовой и жидкой фаз углеводородной системы и объема внедряющейся в залежь воды. Для зоны I, где r_v>r_k (рис. 2), выражения для вычисления текущего дебита скважины напишем в следующем виде (31), (32):

$q_1=-\frac{d}{dt}\left\{\left[\frac{s_1}{B\left(p_1\right)}+\left(1-s_1\right)\frac{p\beta c\left(p_1\right)}{z\left(p_1\right)p_{at}}\right]{\omega }_1\left(p_1,t\right)\right\},$ (31)

$q_{g1}=-\frac{d}{dt}\left\{\left[\begin{array}{l}\frac{\left(1-s_1\right)p_1\beta }{z\left(p_1\right)p_{at}}\times \\ \times \left[1-c\left(p_1\right)\overline{\gamma }\left(p_1\right)\right]+\\ +\frac{s_{1}S\left(p_1\right)}{B\left(p_1\right)}\end{array}\right]{\omega }_1\left(p_1,t\right)\right\},$ (32)

где p₁ и s₁ – средние давление и насыщенность пор жидкостью в зоне I, соответственно; ${\omega }_1={\Omega }_{10}{\overline{m}}_1-\frac{1}{1-{\rho }_{ost}}{\overline{m}}_1{V}_v,$ V_v – объем вторгшейся в залежь воды и текущее положение фронта воды r_v определяются с учетом (10).

А в зоне II (рис. 3), где средневзвешенное давление и насыщенность пор жидкостью соответственно p и s, уравнения материального баланса имеют следующий вид [41] (33), (34):

$q-q_1=-\frac{d}{dt}\left\{\left[\frac{s}{B\left(p\right)}+\left(1-s\right)\frac{p\beta c\left(p\right)}{z\left(p\right)p_{at}}\right]\Omega \left(p,t\right)\right\},$ (33)

$q_g-q_{g1}=-\frac{d}{dt}\left\{\left[\begin{array}{l}\frac{\left(1-s\right)p\beta }{z\left(p\right)p_{at}}\times \\ \times \left[1-c\left(p\right)\overline{\gamma }\left(p\right)\right]+\\ +\frac{sS\left(p\right)}{B\left(p\right)}\end{array}\right]\Omega \left(p,t\right)\right\}.$ (34)

Здесь W – текущий объем пор, насыщенных углеводородами; q, q_g – дебит жидкости и газа скважины; q₁, q_g1 – расход жидкости и газа, перeтекших из зоны II в зону I через границы , когда r_v<r_k уравнения материального баланса выписываются в виде (35), (36):

$q=-\frac{d}{dt}\left\{\left[\frac{s}{B\left(p\right)}+\left(1-s\right)\frac{p\beta c\left(p\right)}{z\left(p\right)p_{at}}\right]\omega \left(p,t\right)\right\},$ (35)

$q_g=-\frac{d}{dt}\left\{\left[\begin{array}{l}\frac{\left(1-s\right)p\beta }{z\left(p\right)p_{at}}\times \\ \times \left[1-c\left(p\right)\overline{\gamma }\left(p\right)\right]+\\ +\frac{sS\left(p\right)}{B\left(p\right)}\end{array}\right]\omega \left(p,t\right)\right\},$

где

$\omega ={\Omega }_0\overline{m}-\frac{1}{1-{\rho }_{ost}}\overline{m}{V}_v.$ (36)

Из (31)–(36) можно получить уравнения, описывающие изменения средневзвешенных пластовых давлений и насыщенностей во времени для периодов вытеснения r_v>r_k и r_v>r_k в следующем виде [42].

Когда r_v>r_k, получим (37), (38):

$\frac{dp}{dt}=\frac{-\frac{q-q_1}{\Omega }\left({\alpha }_4+G{\alpha }_2\right)}{\left({\alpha }_5+{\alpha }_6\right){\alpha }_4+\left({\alpha }_7+{\alpha }_8\right){\alpha }_2},$ (37)

$\frac{ds}{dt}=\frac{-\frac{q-q_1}{\Omega }G-\left({\alpha }_7+{\alpha }_8\right)\frac{dp}{dt}}{{\alpha }_4},$ (38)

$\frac{d\Omega }{dt}={{\overline{m}}^{\text{'}}}_p\frac{dp}{dt},$ если учитывать, что породы коллектора подвергаются нелинейно-эластической деформации, то получим (39)–(44):

$\overline{m}=\exp \left[a_m\left(p-p_0\right)\right],$

то

${{\overline{m}}^{\text{'}}}_p=a_m\text{exp}\left[a_{m}1\right(p-p_0\left)\right].$ (39)

$\frac{d{p}_1}{dt}=\frac{-\frac{q_1}{{\omega }_1}\left({\alpha }_4+G_1{\alpha }_2\right)+\left({\alpha }_2{\alpha }_3+{\alpha }_1{\alpha }_4\right)\frac{q_v}{{\omega }_1}{\overline{m}}_1}{\left(\begin{array}{l}\left({\alpha }_5+{\alpha }_6\right){\alpha }_4+\left({\alpha }_7+{\alpha }_8\right){\alpha }_2+\\ +\left({\alpha }_2{\alpha }_3+{\alpha }_1{\alpha }_4\right)\frac{\left(1-{\rho }_{ost}\right)}{{\omega }_1}{{\overline{m}}^{\text{'}}}_{p1}\left({\Omega }_{01}-\frac{V_v}{1-{\rho }_{ost}}\right)\end{array}\right)},$ (40)

$\frac{d{s}_1}{dt}=\frac{\left(\begin{array}{l}-\frac{q_1{G}_1}{{\omega }_1}-\left({\alpha }_7+{\alpha }_8\right)\frac{d{p}_1}{dt}-{\frac{{\alpha }_3}{\omega }}_1\left(1-{\rho }_{ost}\right)\times \\ \times \left[\frac{d{p}_1}{dt}{{\overline{m}}^{\text{'}}}_{p1}{\left({\Omega }_{01}-\frac{V_v}{1-{\rho }_{ost}}\right)}_1-\frac{q_v}{1-{\rho }_{ost}}{{\overline{m}}^{\text{'}}}_1\right]\end{array}\right)}{{\alpha }_4},$ (41)

$\frac{d{\Omega }_1}{dt}=\left({\Omega }_{01}-\frac{V_v}{1-{\rho }_{ost}}\right)\frac{d{\overline{m}}_1}{dt}-\frac{q_v}{1-{\rho }_{ost}}{\overline{m}}_1,$

где

$\begin{array}{c}{\overline{m}}_1=\exp \left[a_{m1}\left(p_1-p_0\right)\right],\\ \frac{d{\overline{m}}_1}{dt}=a_{m1}\exp \left[a_{m1}\left(p_1-p_0\right)\right]\frac{d{p}_1}{dt},\end{array}$  (42)

и, когда r_v<r_k:

$\frac{dp}{dt}=\frac{-\frac{q}{\omega }\left({\alpha }_4+G{\alpha }_2\right)+\left({\alpha }_2{\alpha }_3+{\alpha }_1{\alpha }_4\right)\frac{q_v}{\omega }\overline{m}}{\left(\begin{array}{l}\left({\alpha }_5+{\alpha }_6\right){\alpha }_4+\left({\alpha }_7+{\alpha }_8\right){\alpha }_2+\\ +\left({\alpha }_1{\alpha }_4+{\alpha }_2{\alpha }_3\right)\frac{1-{\rho }_{ost}}{\omega }{{\overline{m}}^{\text{'}}}_p\left({\Omega }_0-\frac{V_v}{1-{\rho }_{ost}}\right)\end{array}\right)},$ (43)

$\frac{ds}{dt}=-\frac{\left(\begin{array}{l}\frac{qG}{\omega }+\left({\alpha }_7+{\alpha }_8\right)\frac{dp}{dt}+\frac{\left(1-{\rho }_{ost}\right)}{\omega }\times \\ \times {\alpha }_3\left[{{\overline{m}}^{\text{'}\text{'}}}_p\frac{dp}{dt}\left({\Omega }_0-\frac{V_v}{1-{\rho }_{ost}}\right)-\frac{q_v}{1-{\rho }_{ost}}\overline{m}\right]\end{array}\right)}{{\alpha }_4},$ (44)

где газовый фактор для внутренней и внешней зон при соответствующих давлениях (p, p₁) и насыщенностях жидкой фазы (s, s₁) определяется следующим выражением [43]:

$G=\frac{\frac{\overline{\mu }\left(p\right)B\left(p\right)p\beta }{z\left(p\right)p_{at}}\left[1-c\left(p\right)\overline{\gamma }\left(p\right)\right]+\frac{S\left(p\right)}{\psi \left(s\right)}}{\frac{1}{\psi \left(s\right)}+\frac{\overline{\mu }\left(p\right)B\left(p\right)p\beta c\left(p\right)}{z\left(p\right)p_{at}}},$ (45)

$\overline{m}=\frac{m}{m_0}$  – отношение текущей пористости внутренней зоны к его первоначальному значению, ${\overline{m}}_1=\frac{m_1}{m_{10}}$  – отношение текущей пористости внешней зоны к его первоначальному значению;

$\begin{array}{c}{\alpha }_1=s\frac{1}{B\left(p\right)}+\left(1-s\right)\frac{p\beta c\left(p\right)}{z\left(p\right)p_{at}},{\alpha }_2=\frac{p\beta c\left(p\right)}{z\left(p\right)p_{at}}-\frac{1}{B\left(p\right)},\\ {\alpha }_3=s\frac{S\left(p\right)}{B\left(p\right)}+\left(1-s\right)\frac{p\beta }{z\left(p\right)p_{at}}[1-c\left(p\right)\overline{\gamma }\left(p\right),\\ {\alpha }_4=\frac{S\left(p\right)}{B\left(p\right)}-\frac{p\beta }{z\left(p\right)p_{at}}\left[1-c\left(p\right)\overline{\gamma }\left(p\right)\right],\\ {\alpha }_5=\left(1-s\right){\left[\frac{p\beta c\left(p\right)}{z\left(p\right)p_{at}}\right]}^{\text{'}},{\alpha }_6=s{\left[\frac{1}{B\left(p\right)}\right]}^{\text{'}},{\alpha }_7=s{\left[\frac{S\left(p\right)}{B\left(p\right)}\right]}^{\text{'}},\\ {\alpha }_8=\left(1-s\right){\left\{\frac{p\beta }{z\left(p\right)p_{at}}\left[1-c\left(p\right)\overline{\gamma }\left(p\right)\right]\right\}}^{\text{'}};\end{array}$

«′» – означает производную по p.

Отметим, что параметры a₁–a₈ вычисляются при соответствующих давлениях и насыщенностях в зависимости от зоны [44–46].

Системы уравнений (37)–(40), (42)–(44) совместно с (24), (25) и (27) с учетом (29) решаются одним из численных методов и позволяют прогнозировать основные показатели процесса вытеснения углеводородной системы, закачиваемой в пласт водой, при условии заданной депрессии или заданного забойного давления [47–49].

Таким образом, получены системы дифференциальных уравнений (37), (38) и (40)–(42) при r_v>r_k, а при r_v<r_k – (43), (44), решения которых при известном дебите скважины q_g позволяют определять среднепластовые давления и конденсатонасыщенности в соответствующих частях коллектора в любой момент времени [50–52].

Результаты исследования и их обсуждение

С учетом решения задачи получен алгоритм для прогнозирования разработки основных показателей процесса вытеснения углеводородной системы к скважине, когда вблизи скважины (во внутренней зоне) и в отдаленной части залежи (т. е. во внешней зоне) пласт имеет разные коллекторско-емкостные и реологические характеристики [53–55].

Вышеизложенный подход позволяет определять основные показатели разработки газоконденсатной залежи при технологических режимах с учетом разных значений в проницаемости и характере деформаций призабойной зоны и отдаленной от забоя части пласта-коллектора [48]. При этом можно использовать нижеприведенный алгоритм.

1. Вводятся начальные значения переменных t=0, p_s=p₀, p_k=p₀, p₁=p₀, s=s₀, s₁=s₀, m=m₀, k₁=k₀₁, k₂=k₀₂, q₁=0 и исходные данные.
2. Вычисляются начальные значения газоконденсатного фактора G₀ по формуле

$G_0=\frac{1-c\left(p_0\right)\overline{\gamma }\left(p_0\right)}{c\left(p_0\right)},$

газонасыщенный объем залежи, запасы газа и конденсата (при s₀=0):

$V_{zap}=\pi R_k^{2}h{m}_0\left[\frac{1-c\left(p_0\right)\overline{\gamma }\left(p_0\right)}{z\left(p_0\right)p_{atm}}\right],V_{kzap}=\frac{V_{zap}}{G_0}.$

3. Определяется текущее положение ГВК (газо-водяного контакта):

3.1) по (10) с учетом (11) и (12) вычисляется текущее значение q_v;

3.2) с помощью выражений (13) определяется объем вторгшейся в залежь воды за промежуток времени Dt и положение ГВК r_v.

4. Если рассматривается случай заданного темпа отбора газа (n процентов в год от начальных балансовых запасов), дебит определяется выражением $q_g=\frac{V_{zap}n}{100}.$
5. Если рассматривается случай депрессии, то задается значение депрессии Dp, иначе переход к шагу «7».
6. Забойное давление вычисляется выражением p_s=p₁–Dp.
7. Рассчитывается дебит скважины.

Если r_v>r_k то:

7.1) вычисляются значения f_g(p, s) по (16*) для давлений p_s, p_v и a, b по (30);

7.2) определяется фиктивная депрессия H_v–H_k по выражению (29);

7.3) вычисляется текущее значение дебита газа по (27) и конденсата $q_k=\frac{q_g}{G};$

иначе:

7.1) вычисляются значения f_g(p, s) по (14) для давлений p_s, p_k, p_v и a, b по (30);

7.2) определяется фиктивная депрессия H_v–H_k, H_k–H_s и H_v–H_s по (29);

7.3) вычисляется текущее значение дебита газа по (24) с учетом (25) и конденсата $q_k=\frac{q_g}{G}.$

8. Вычисляются текущие значения конденсатонасыщенности и среднепластового давления для времени t+Dt.

Если r_v>r_k:

Текущие значения p, s во внешней зоне определяются решением системы дифференциальных уравнений (40)–(42), а во внутренней зоне – системы уравнений (37)–(39).

Иначе:

Текущие значения  во внутренней зоне определяются численным решением системы уравнений (43), (44).

9. Вычисляется текущее значение газоконденсатного фактора G по (45).
10. Определяются текущие значения накопленного отбора газа и конденсата и, следовательно, их коэффициенты извлечения:

$K_g=\frac{{\displaystyle {\sum }_{t=0}^t{q}_g\Delta t}}{V_{zap}}\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}и\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{K}_k=\frac{{\displaystyle {\sum }_{t=0}^t{q}_k\Delta t}}{V_{kzap}}.$

11. Проверяется значение пластового давления, если оно больше заданного его значения как конечного, переход к пункту «3», иначе переходим к шагу «12».
12. Вывод результатов и конец.

Заключение

Полученное выше решение позволяет прогнозировать основные показатели разработки газоконденсатной залежи, представленной упругими зонально-неоднородными коллекторами. При этом учитываются реальные PVT свойства двухфазной углеводородной системы и реология пород-коллекторов. Предложенный алгоритм позволяет моделировать практически любой технологический режим закачки и скважины. Так, можно воспроизвести режим заданного темпа закачки воды и заданного давления на контуре заводнения. Следует отметить, что в случае, если принять темп закачки равным нулю, получается модель истощения. Относительно режима скважины отметим, что можно смоделировать случай, когда задается депрессия, а также можно моделировать режим заданного забойного давления и заданного отбора. В последнем случае не приходится вычислять дебит.

Авторлар туралы

Marat Khabibullin

Ufa State Petroleum Technical University

Хат алмасуға жауапты Автор.
Email: m-hab@mail.ru
ORCID iD: 0000-0003-2565-0088

Cand. Sc., Associate Professor

Ресей, Oktyabrsky

Rustem Suleimanov

Ufa State Petroleum Technical University

Email: burmash2017@mail.ru
ORCID iD: 0000-0003-2878-9253

Cand. Sc., Associate Professor

Ресей, Oktyabrsky

Razifa Stepanova

Ufa State Petroleum Technical University

Email: razifa0210@yandex.ru

Cand. Sc., Associate Professor

Ресей, Oktyabrsky

Alina Gizzatullina

Ufa State Petroleum Technical University

Email: alina.gizzatullina87@mail.ru
ORCID iD: 0000-0002-8719-9942

Cand. Sc., Associate Professor

Ресей, Oktyabrsky

Arsen Khabibullin

Ufa State Petroleum Technical University

Email: zx9ooqpv@gmail.com

Student

Ресей, Oktyabrsky

Әдебиет тізімі

Қосымша файлдар