Ключевые слова:
математическое моделирование, газоконденсат, фазовые переходы, фильтрация в пористой среде.
Keywords:
mathematic modeling, gas condensate, changes of states, filtration in porous medium.
УДК 622.276.03:532
В.В. Качалов, И.Л. Майков, ДА Молчанов, В.М. Торчинский
Газоконденсатная залежь как колебательная система осцилляторного типа
Смесь углеводородов газоконденсатных месторождений представляет собой сложную смесь метана и высших производных метанового ряда с высоким содержанием метана. Поведение газоконденсатной системы как в пласте на значительном уда -лении от скважины, так и вблизи нее определяется большим количеством факторов: исходным составом смеси, температурой, пластовым давлением, градиентом давления, свойствами породы коллектора. На процессы фильтрации влияют вязкость фаз, поверхностное натяжение, краевые углы, структура коллектора. Фазовая диаграмма смеси углеводородов содержит так называемую ретроградную область, в которой при снижении давления вплоть до давления максимальной конденсации возможно образование ретроградной жидкости, испаряющейся с дальнейшим уменьшением давления [1, 2]. Следует отметить, что наряду с нетипичной фазовой диаграммой для газоконденсатных смесей характерны дополнительные особенности в поведении функций фазовых проницаемостей, которые проявляются при фильтрации смеси углеводородов через пористую среду. Значения функций относительных фазовых проницаемостей и для газовой, и для жидкой фаз могут иметь нулевые значения.
Неустойчивость работы газоконденсатных месторождений иногда отмечалась на скважинах и в некоторых экспериментальных исследованиях [2-4]. Однако природа этих неустойчивостей остается неясной [3].
Целью настоящей статьи является теоретическое исследование особенностей фильтрации смеси углеводородов в пористой среде на примере бинарной смеси метан - н-бутан.
Математическая модель
Рассматривается одномерная нестационарная фильтрация двухфазной двухкомпонентной системы (метан - н-бутан) в изотермических условиях в предположении равенства давления в фазах при условии фазового равновесия (характерные времена фазовых переходов существенно меньше гидродинамических).
Фазовая диаграмма такой смеси [5, 6] показана на рис. 1. Состояние 1 на фазовой диаграмме соответствует входному граничному условию P1, состояние 2 - выходному граничному условию Р2. За счет фильтрации смесь углеводородов переходит из состояния 1 в состояние 2. Обозначим длину, на которой давление изменяется от Р1 до Р2, как L.
Процесс фильтрации при допущении о химической нейтральности компонентов в одномерном случае описывается уравнениями баланса массы для каждого компонента, записанными в дивергентной форме относительно молярных плотностей [3]:
ш-
dt
Е
V J
CijnjSj
дх
Е c»nPi
V 1
= 0,
(1)
где m - пористость; индекс i соответствует компоненту (CH4, C4H10); индекс j соответствует фазе (1 - газ, 2 - жидкость); ctj - мольная доля i-го компонента в j-й фазе двухфазной смеси; n - молярная плотность j-й фазы; sj - объемная доля j-й фазы в смеси; Uj - скорость j-й фазы; t - время; х - координата.
№ 2 (18) / 2014
Актуальные вопросы исследований пластовых систем месторождений углеводородов
107
Рис. 1. Фазовая диаграмма бинарной смеси метан - н-бутан при температуре 310 K: c1(P) - кривая равновесных концентраций в газе; c2(P) - кривая равновесных концентраций в жидкости
Уравнения сохранения импульса записываются в приближении закона Дарси:
Uj = —к
fj (Sj ) dP Ц j dx ’
(2)
где к - абсолютная проницаемость; у. - коэффициент динамической вязкости j-й фазы; P - давление;f - функция относительной фазовой проницаемости j-й фазы.
Запишем уравнения состояния для газовой и жидкой фаз, выраженные через коэффициенты сжимаемости [2, 6]:
P = njZjRT, (3)
где z. - коэффициент сжимаемости j-й фазы; R - газовая постоянная; T - температура смеси.
Равновесные концентрации компонент в обеих фазах зависят от давления и удовлетворяют условию
X Cj( P) = 1. (4)
i
Объемные доли фаз связаны соотношением
X ^ = 1. (5)
J
Система уравнений дополняется граничными и начальными условиями.
Численный метод
При численном решении системы уравнений необходимо учитывать наличие особой точки, где существуют разрывы некоторых функций. Введем функцию газонасыщенности s (тогда влагонасыщенность запишется как (1 - s)), а также следующие обозначения: c1 (P) = cCH д -мольная доля метана в газе; c2(P) = cCH 2 -мольная доля метана в жидкости (тогда мольная доля бутана в газе выразится как (1 - с1), а в жидкости - (1 - с2)).
С учетом принятых обозначений систему (1)-(2) запишем в виде
m—(c1n1s + c2 n2 (1 — s)) —
d f, ( fi(s) + f2(s) 1 ЭР', 0
— I k\ cini------+ c2n2 ----- I— 1 = 0,
dx I ^ Ц1 Ц2 J dx
m—((1 — c1)n1s + (1 — c2 )n2 (1 — s)) — dt
d (k ((1 ) f1(s) , (1 ) f2(s) 1 dP
— Г— I к I (1 — c1 ) n1--+ (1 — c2 )n2 ---- I —
dx I ^ Ц1 Ц 2 J dx
0. (6)
Введем переменные
Y, = c1n1s + c2n2(1 - s),
Y2 = (1 - cl)nls + (1 - c2)n2(1 - s). (7)
№ 2 (18) / 2014
108
Научно-технический сборник • ВЕСТИ ГАЗОВОЙ НАУКИ
Здесь новые переменные представляют собой количество молей метана и бутана в газовой и жидкой фазах, причем введенные таким образом функции непрерывны.
С учетом уравнений (7) и (3) систему уравнений (6) можно записать как
ia If < т+c т 1 р эр 1
dt mRT Эх lljZj _2 z2 J dx J
dlL
dt
k Э mRT dx
(1 - ci)
fi(s)
_ z1
+ (1 - Cl)
f2(s) _ 2 z2
(8)
Система (8) описывает изменение количества молей метана и бутана во времени и пространстве.
Складывая уравнения системы (6) и подставляя выражения для плотности из уравнения (3), получим уравнение для изменения по-рового давления:
ЭРВ22 к Э (ЛМ) , fi(s) 1
Э/ m Эх ^ 1 _1 г _2 Z2 у
Г. 1 - ,v ^
где D22 = —+
1Y Z2 у
р ЭР
Эх
Л
(9)
Запишем систему уравнений (8)-(9) в безразмерном виде:
^=± гг ст+ст) р ip
Эт ЭХ ^ _zj _2z2 ) dX j
^-± Ц и - cj) f«+и - с.) 1 р ±
Эт ЭХ И _Zj z2 1 ЭХ I
ЭрРгг = Э ГГ Ms) + f.js) | Эр
Эт ЭХ |l pjZj p.z. J ЭХ I’
(10)
Y
Y,
= P . p _ P.
где y = —; y2 = —; Y0 = -Y-; p = —; X = -;
Y
t
т = —; то _
Y
^Pi( P) L\ _
RT
P
L
_i _■
_1 -77
_2
_2
т0 ' kPi _i(Pi) ' _l(Pl)
Функция газонасыщенности определяется через значения переменных ^ и у2:
С - с2(р)
s = -
с1(Р) - c2(р) где мольная доля метана в смеси
(11)
C = У1 . (12)
У1 + У2
Правые части системы уравнений (10) с учетом (11) и (12) являются функциями p, у1 и у2. Система уравнений (10)-(12) замыкается
начальными и граничными условиями для переменных y1, у2 и p. При численном решении используется полностью неявная схема для всех указанных переменных.
Приближенный анализ
Для решения системы уравнений (8)-(11) необходимо дополнительно иметь функциональные зависимости c1(p), c2(p), z1(p) и z2(p), которые вычисляются из условий фазового равновесия [5, 6]. Вид кривых Cj, c2 для смеси метан -н-бутан отражен на рис. 1. Функциональные зависимости вязкостей ц1(р), ц2(р) вычисляются согласно алгоритму, представленному в работе [6]. Функции c2(P), z1(P), z2(P), ц1(Р) и ц2(Р) являются непрерывными и монотонными на отрезке [P2; P1], а функция c1(P) - унимодальной. Функциональные зависимостиf(s) иf2(s) имеют вид [6, 7]:
Ms)
fiis)
0, если 0 < s < s1k;
[(s - sik ) 1 (1- sik )f'5 ' если sik < s < i;
[(s2 k- s) 1 s2k ]3'5' если 0 <s < s2k;
0, если s2k < s < 1,
(13)
где s1k и s2k - параметры. В работе принимались значения s1k = 0,35 и s2k = 0,8. Основная особенность зависимостей (13) состоит в том, что существуют области, в которых каждая из функций становится равной 0, т.е. скорость ее движения равна 0. Физически это означает возникновение областей, в которых фаза не является непрерывной, а присутствует в виде пузырей или капель, которые неподвижны.
Рассмотрим упрощенную постановку задачи. Предположим, что вязкости фаз постоянны и
1
_2( Ф
<< 1,
(14)
коэффициенты сжимаемости равны 1. Тогда систему уравнений для вычисления у1 и p (см. (10)) можно представить в виде (в предположении Ар = 1 - р2 << 1):
^ = _± Эт ЭХ
s) + а1с2 /2( s)) ЭХ
Эр = _ Э Эт = _ ЭХ
(f( s) + а/>(s))
Э_
ЭХ
(15)
№ 2 (18) / 2014
Актуальные вопросы исследований пластовых систем месторождений углеводородов
109
Заменим зависимости относительных фазовых проницаемостей (13) на модельные функции 0:
fi(s) = kfi(s - su),
f2(s) = £2(1 -0(s - s2t)), (16)
Решение задачи (20) дает значение давления в точке Xp = AX:
1 - Xp + aXp (1 -Ap)
= 1 - Xp + aXp
k2
где a = a, —.
k
(21)
где k1 < 1; k2 < 1.
Пусть значение газонасыщености находится внутри интервала (s1k, s2k), т.е. f(s) = k1 и f2(s) = k2. Рассмотрим режим течения, который установился через достаточно большой начальный промежуток времени. Систему уравнений (15) можно записать в виде
Таким образом, возникает зона, в которой давление близко к входному граничному, но
s < s1k.
Пусть следующая точка, в которой s достигла величины s1 k, находится на расстоянии Xp = AX, = 2 AX. Тогда установившееся давление определяется из решения задачи
М = - ( др_
дт I дХ
d2 р
~dX
= 0.
К
—L + a,k2
dp
dc2
dp
(17)
Таким образом, процесс фильтрации удо-влетворяет системе уравнений (17): в каждый момент времени рассматриваем установившийся режим по давлению и нестационарный процесс фильтрации. Решение для давления имеет вид
p = 1 - ApX.
(18)
Согласно первому уравнению системы (17) и с учетом (14):
м
дт
ki p
х Y х < о
dX J dp
(19)
Так как max
max
дт
dcj
dp
достигается при X ^ 0, то и
достигается при X ^ 0. Уменьшение
у, приводит к уменьшению C и s (согласно (12) и (11)). Пусть за время Ат газонасыщенность s в точке AX (AX << 1) уменьшилась до s1k, т.е. fj(s) = 0. Установившееся давление удовлетворяет следующим уравнениям:
d 2 Pi
dX2
d 2 Pii
dX2
0, pi( X = 0) = 1,
0, Pii( X = 1) = 1 -Ap,
Pi( X = Xp) = Pii( X = Xp),
k dpHX = Xp) = a,k2 dpL(X = Xp), 1 dX p 12 dXy p
где Xp = AX.
(20)
d 2 Pi dX2
0, Pi( X = 0) = 1,
—p- = 0, pII( X = 1) = 1 -Ap,
dx2
Pi( X = Xp) = Pii( X = Xp),
addp^(X = X ) = dPp^{X = X ),
dX p dX p
(22)
а давление в точке Xp - выражением
a + (1 -Ap)Xp -aXp a + Xp-aXp
(23)
Выражение (23) последовательно дает решение для фронта давления в точках Xp = (n + 1) и AX(n > 1), в которых газонасыщенность достигает величины s1k (т.е. процесс фильтрации является последовательностью стационарных процессов (по давлению), отличающихся положением фронта давления, на котором газонасыщенность достигает величины s1k).
Зависимости (21) и (23) pp от Xp представлены на рис. 2.
Аппроксимируем линейной зависимостью зависимостьpp отXp (рис. 3):
Pp = 1 -вXp,Xp <Xb,
Pp = 1 -Ap + у(1 - Xp), Xp > X4,
X = ap - y
b P-y ,
p,=i-P^Ap. P-y
(24)
Рассмотрим движение по прямой ab. Введем обозначения: у0 - значениеу,, при котором газонасыщенность достигает критического значения s = s1k; у'" - значение у, при т = 0, у'” > у0. Точки на прямой ab соответствуют времени достижения критического значения
№ 2 (18) / 2014
110
Научно-технический сборник • ВЕСТИ ГАЗОВОЙ НАУКИ
Рис. 2. Зависимостьpp от Xp при а = 0,01 (1) и а = 0,05 (2)
Рис. 3. Аппроксимационная зависимость pp от Xp
газонасыщенности s1k. Новое распределение давления ab1c достигается в два шага:
1) согласно уравнению (19), у1 за время т достигает значения yj0. Изменение у1 определяется по формуле
2) происходит изменение давления на величину Дд (Xb^) = p - p (достижение точки b1). Полное изменение у1 определяется по формуле (при этом уже s < s1k, соответственно,
yi < у0)
ДУ1X) = pKxl |р (Ра, К < о, (25)
1 dp
dp
где Xi - градиент давления на прямой ac;
ЭХ
ДУ X) = pKxl (Pa, )T1 + с1ДРь X) < ° (26)
Пусть в момент времени т2 распределение давления соответствует кривой ab2c.
№ 2 (18) / 2014
Актуальные вопросы исследований пластовых систем месторождений углеводородов
111
Изменение у, за время т2 определяется по формуле (согласно (19) и (26))
) = frxf d-iраг ^ +
+ С1 ЛРе2 (Х2 ) + Р^1%2 (ре2 )(Т2 - X! ), (27)
где ЛРе2 (ХЬ2 ) = Pa2 - Pes Ъ - ГРаДиент Давления
др „ ,
дХ на прямой b,0.
Первые два члена отвечают за изменение у, за время т,, третий член - за изменение у, за время т2 - т,, при этом газонасыщенность достигает критического значения s = s,k (соответственно, у, = у0).
Полное изменение у, определяется по формуле (с учетом изменения давления)
— 2 дс,
А» (Xb2 ) = PKh -Х-(Ра2 К +
др
дс,
+ PkX Эр (Pe2 )(Т2 -Т1) +
+ с,(Лрр2(X) + Лр(,X )),
(28)
гДе APb2( Xb2) = Pb2 - Ppe2.
Таким образом, можно последовательно получить время достижения любой точки распределения давления pp (например, распределение abno, см. рис. 3).
Система для определения Ayl(Xb ) упро-
дс1 , ч
щается, если предположить, что —L (pn) = ф =
dp n
= const и = X2, причем |у| < |Х| < |Др|. Тогда для момента времени тп можно записать следующее уравнение:
Ay (X„) = pklX2^in + cxApn (Хп), (29)
где Хп = (п + ^ДХ, п > ! и Ар (X) - {X^-р)ХsX•
1(1 -X)(1-Ap),X, >X,. (30)
Уравнение (29) с учетом (П), (!2) можно записать в виде
Тп =~ГРТ~ I Asn (С ( Pn ) - С2 ( Pn )) - 4 C1 ( Pn )APn (Xn ) | , (3 , )
где Дsn = s,k - sn. В формулах (29) и (3!) 1 = у
при Хп > Xb.
Физически точка b разделяет два режима течения. Поток в точках b 2 прямой ab лимитируется условиями выше по течению (проводимость выше по потоку определяется свойствами жидкой фазы, ниже - свойствами газовой фазы). Поток в точках о2 прямой bo лимитируется условиями ниже по потоку (и хотя проводимость ниже по потоку определяется свойствами газовой фазы, величина потока ограничивается малым градиентом давления на прямой bo).
При достижении величины газонасыщенности s < s№ согласно (26) или (28), начинается процесс испарения (увеличение у,):
— = ак, />PV > 0, Эх
Эс2
где ^ = ~TLiPn). dp
(32)
При движении по прямой bo (точки о, и o2, см. рис. 3) в точках прямой ab возникают дополнительные члены, связанные с положительным скачком давления (например, распределению давления ao, o на рис. 3 соответствует Apn = ph - pb > 0), т.е. общее изменение у, на этапе испарения описывается формулой, аналогичной (29):
= akрР V + с\Лрп > 0. Эх
(33)
При этом согласно (34) происходит увеличение у, и s > s1k. При достижении правой границы профиль давления восстанавливается.
Рассмотрим возможные режимы фильтрации. Движение в системе осуществляется за счет заданного перепада давления, конкретные значения давления определяют участок фазовой диаграммы и все остальные свойства колебательной системы. За счет внутренних свойств системы происходит уменьшение расхода: выпадение конденсата в ретроградной области, уменьшение значения фазовой проницаемости газа и проводимости выше по потоку (переход с прямой ao на прямую ab) или уменьшение расхода за счет уменьшения градиента давления (переход с прямой ao на прямую bo). В рассмотренном выше случае газовая фаза останавливается. Хотя плотность жидкой фазы выше газовой, скорость ее движения существенно ниже скорости газовой фазы (см. (,4)). Поэтому энергия системы уменьшается (стадия возврата энергии во внешний источник) и достигает некоторого минимального
№ 2 (18) / 2014
112
Научно-технический сборник • ВЕСТИ ГАЗОВОЙ НАУКИ
значения, при этом AE™ - максимальное количество отданной энергии за период. Затем начинается увеличение расхода (испарение выпавшего конденсата, увеличение фазовой проницаемости газа и увеличение проводимости выше по потоку). Энергия системы начинает возрастать за счет вовлечения в движение остановившейся газовой фазы (стадия отбора энергии из внешнего источника) и достигает некоторого максимального значения, при этом AE““ - максимальное количество подведенной энергии из внешней системы за период. Если
АСГ -AEim“ > 0, (34)
то в течение каждого периода колебаний энергия у системы отбирается, и происходят затухающие колебания. Если наоборот:
ACT -AET < 0, (35)
то энергия системы возрастает, и колебания возрастают. Таким образом, система представляет автоколебательную систему осцилляторного типа [8]. Здесь имеется источник питания системы энергией, причем подвод энергии происходит не произвольно, а осуществляется при помощи механизма управления, приводимого в действие самой системой. Механизм управления действует как обратная связь между колебательной системой и источником энергии,
Список литературы
1. Вяхирев Р.И. Разработка и эксплуатация газовых месторождений / Р.И. Вяхирев,
A. И. Гриценко, Р.М. Тер-Саркисов. - М.: Недра, 2002. - 880 с.
2. Зайченко В.М. Моделирование процессов фильтрации углеводородов в газоконденсатном пласте / В.М. Зайченко, И. Л. Майков,
B. М. Торчинский и др. // Теплофизика высоких температур. - 2009. - Т. 47. - № 5. - С. 701-706.
3. Николаевский В. Н. Геомеханика и флюидодинамика / В.Н. Николаевский. - М.: Недра, 1996. - 447 с.
4. Митлин В. С. Подземная гидромеханика сложных углеводородных смесей /
B. С. Митлин // ВИНИТИ. - 1991. - Т. 4. -
C. 154-222.
обеспечивая подвод энергии в нужный момент периода колебаний.
Проведенные исследования фильтрации смеси углеводородов показали, что возникающие колебания описываются колебательной системой осцилляторного типа. Необходимым условием возникновения колебаний является нахождение газоконденсатной системы в ретроградной области на фазовой диаграмме, когда при понижении давления смеси происходит выпадение жидкости. Свойства колебательной системы однозначно определяются граничными условиями для давления (выделяется область по давлению на фазовой диаграмме смеси) и мольной долей легкого углеводорода на входе (вычисляется конкретная точка на фазовой диаграмме смеси). Наличие областей с нулевой фазовой проницаемостью является основой осуществления обратной связи колебательной системы. В зависимости от свойств колебательной системы возможны как затухающие, так и незатухающие колебания -возникновение автоколебаний. Таким образом, колебательная система осцилляторного типа в зависимости от условий может реализовывать различные режимы фильтрации как общей системы, так и ее отдельных элементов. Так, в промысловых условиях реализуется режим образования зон неподвижного ретроградного газового конденсата или защемленной газовой фазы, обогащенной метаном.
5. Баталин О.Ю. Фазовые равновесия в системах природных углеводородов / О.Ю. Баталин,
A. И. Брусиловский, М.Ю. Захаров. - М.: Недра, 1992. - 272 с.
6. Директор Л.Б. Одномерная нестационарная модель двухфазной фильтрации газоконденсатной смеси / Л.Б. Директор,
B. В. Качалов, И. Л. Майков и др.; препринт // ОИВТ РАН. - М., 2000 - № 2 (441). - 46 с.
7. Вафина Н.Г. К определению фазовых проницаемостей в случае многокомпонентной фильтрации с фазовыми переходами /
Н.Г. Вафина, С.Н. Закиров, П.А. Юфин // Известия вузов. - 1988. - № 9. - С. 59-63. -(Серия «Нефть и газ»).
8. Магнус К. Колебания: введение в исследование колебательных систем / К. Магнус; пер. с нем. -М.: Мир, 1982 - 304 с.
№ 2 (18) / 2014