Моделирование процессов фильтрации в коллекторах с переменной пористостью Текст научной статьи по специальности «Математика»

УДК 519.688

А. В. Шевченко, И. В. Цыбулин, Ю. И. Скалько

Московский физико-технический институт

Моделирование процессов фильтрации в коллекторах с переменной пористостью

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

Ключевые слова: многофазная многокомпонентная фильтрация, переменная пористость, кероген.

A.V. Shevchenko, I. V. Tsybulin, Y.I. Skalko

Moscow Institute of Physics and Technology (State University)

Modeling of the filtration process in reservoirs with

variable porosity

The accountment of chemical reactions that leads to a significant change in porosity is crucial for adequate modeling of thermal EOR methods on shale and Bazhenov Formation reservoirs. In this paper, we construct a mathematical model and a computation algorithm designed for simulation of these problems. The simulator is compared with its commercial analogs. Examples of calculations with variable porosity are given.

Key words: multiphase multicomponent filtration, variable porosity, kerogen.

1. Введение

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

Целью данной работы является разработка математической модели, способной учитывать процессы изменения пористости при реализации тепловых методов воздействия на пласт. Разработанная модель реализована в вычислительном комплексе. Приводится сравнение расчётов с коммерческим симулятором Schlumberger Eclipse для верификации, после чего приведены расчёты с учётом изменения пористости.

2. Математическая модель

Для упрощения математической модели принимаются следующие положения.

• Задача рассматривается в двумерной постановке в плоскости (х, z) (вертикальное и латеральное направления).

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

• Возможны фазовые превращения (испарения лёгких нефтей и Н2О).

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

• Все процессы являются локально квазистатическими (или квазиравновесными) [1]. Это означает, что в каждой расчётной ячейке термодинамическое равновесие устанавливается очень быстро.

Разработанный симулятор предполагает произвольную структуру и физическо-химические свойства резервуара в вертикальном и латеральном направлениях, в том числе переменную пористость, проницаемость, начальный компонентный состав, давление и температуру.

2.1. Уравнение баланса количества вещества

Уравнение баланса количества вещества для компонента г может быть записано в виде [2]:

дИ

где

дМг +ё1удг = (1)

N1 — ^^ Хг>апава — молярная плотность компонента г,

а

• ^ = хг,аПа^а — поток молярной плотности компонента г,

а

• Si — источники/стоки г-го компонента (состоят из источников/стоков химического и скважинного типов),

• %г,а — молярная концентрация компонента г в фазе а (то есть количество молей компоненты г в фазе а в расчёте на один моль фазы а). Суммирование по всем компонентам даст единицу: ^.¡^ — 1,

• па — молярная плотность фазы а,

• 9а — объёмная доля фазы а.

Скорости фильтрации подвижных фаз описываются законом Дарси [3]:

W« — -^(ур - ра%), (2)

где

• ка — относительная фазовая проницаемость фазы а,

• Ц-а — динамическая вязкость фазы а,

• ра — массовая плотность фазы а,

• К — тензор абсолютной проницаемости (диагональный).

Уравнения баланса (1) записываются в том числе и для неподвижных компонентов (скелет, кероген, кокс). В этом случае, поскольку для неподвижной фазы 5 скорость фильтрации полагается равной нулю: = 0, молярная плотность неподвижных компонентов может изменяться только за счёт химических реакций.

Пористость является производной величиной и вычисляется как сумма объёмных долей подвижных фаз:

V = 0Ь + 0С + вш. (3)

Из условия (3) насыщенности подвижных фаз рассчитываются как

0ь Ос

8Ь = —, 8С = -, = -. (4)

V V V

Исходя из (4), насыщенность скелета рассчитывается как вз = 1 — <р = 1 — 8ь — во — вщ.

2.2. Уравнение баланса энергии

Уравнение баланса энергии имеет вид

дЕ

— = С, (5)

где

Е = ^ па0аеа = ^2 Па— р — плотность энергии (см. [4])

а а

J = ^^ паНаШа — каваУТ — плотность потока энергии,

а а

• еа — молярная энергия фазы а,

• На — молярная энтальпия фазы а,

• С — источники/стоки скважинного типа.

Тепловые эффекты фазовых переходов включаются в молярные энергии (и энтальпии) фаз.

2.3. Замыкающие соотношения

Уравнения баланса количества вещества (1)—(2) и уравнение энергии (5) являются основными уравнениями модели (всего п + 1 уравнение, где п — число компонентов). При этом неизвестными переменными модели являются количество молей компонентов и температура (п + 1 неизвестная). Однако для определения величин, входящих в указанные уравнения, требуются дополнительные замыкающие соотношения.

• Поскольку возможны фазовые переходы, а процессы фильтрации значительно медленнее термодинамических процессов, многофазная смесь должна находиться в равновесном состоянии. Это означает, что должны выполняться условия фазового равновесия в рассматриваемой системе. Более подробно с учётом фазового равновесия в вычислительном методе можно ознакомиться в [5].

• Должны быть определены зависимости, используемые в уравнениях баланса: Н% = кг(р,Т) (молярная энтальпия), Ъг^с — Ы,ь = ДН = ДЫ(р,Т) (теплота фазового перехода для испаримых компонентов), ьа = Уа(р,Т,Х1}01) (уравнения состояния фаз), ра = Ра(з>,Т,х^а) (корреляции вязкости), па = па(Хг,а) (относительные фазовые проницаемости).

• Должны быть заданы источники/стоки вещества и энергии скважинного типа (например, можно использовать зависимости, задаваемые формулой Писмана [6]).

• Должны быть заданы источники/стоки химического типа.

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

Разложение керогена

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

где т — характерное время разложения керогена (например, 100 с), Икег,о — начальная концентрация керогена, £(Т) — 1 — а(Т), а а(Т) — часть керогена, разлагающаяся при температуре Т (монотонно возрастающая от 0 до 1 функция), ^,5 — генерационные потенциалы керогена (количество моль нефти и кокса, выделяющихся на один моль разложившегося керогена).

3. Метод решения

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

• р — давление,

• Сг — относительные концентрации компонентов,

• ц — молярная энтальпия смеси.

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

Всего было определено п + 1 уравнений, описывающих поведение системы. При этом основных величин, указанных выше, п + 2. Это является допустимым, поскольку должно выполняться ещё одно соотношение:

Кег ^ 7011 + 50,

^ = — 1(^кег — N^,0 • йТ)),

^ = -№ег — ^Кег,0 • £(Т)), т т

^ = - —N^,0 • т),

.ег

(6)

3.1. Уравнение пьезопроводности

Уравнение пьезопроводности описывает эволюцию давления в пласте. Данное уравнение является следствием из основных уравнений модели, и идеи, используемые в его выводе напрямую используются в вычислительном алгоритме.

Как было сказано ранее, параметры смеси (состав и доли фаз, температура) полностью определяются по основным величинам р, д, ц, то есть

= Иг(р, С1,С2, ...,сь,Г1), г = !,...,П, Е = Е(р,С1,С2,.. .,сь,гО,

где п — количество компонентов. Отсюда

дЩ _ дЫг др дЩ д^ дЩ &п

др дсу дц '

дЕ дЕ др дЕ дс^ дЕ дг/

др дЪ дсу д^ '

(7)

Учтем, что = Ид, Е = Н — р = И-ц — р. Пусть z = (д, с2,..., сь, ц). Таким образом, уравнения (7) можно переписать в виде

дИг _ др д(Ид) д-р д(Ид) дс^ д(Ид) &п дЪ др дсу д^

дЕ _ др д(N71) др д(И^) дс^ д(И^) д^ (8)

д1 др дЪ дсу д^ '

Возьмем такую линейную комбинацию уравнений (8), чтобы из них исключались все производные, кроме др/дЪ. Такая же линейная комбинация строк матрицы д(/дNz)z должна обратиться в ноль, поскольку именно эта матрица стоит перед производными дz/дt.

Утверждается, что такой линейной комбинацией будет совокупность уравнений с коэффициентами вектора дV/дz. Действительно,

у^ дУ д (Яъ) = + у^ дУ = мдУ + д^ у^ дУ ^ = мдУ + ду_

В последнем равенстве д означает показатель однородности функции V(р, z) по аргументам z, то есть

У(р,аь) = аяУ(р, z).

Поскольку У определена только на многообразии ^^ д = 1, ее можно произвольно продолжить в область ^г д = 1. Предлагается продолжить функцию V так, чтобы она была линейно-однородна по z (д = 1):

сг сг

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

Возвращаясь к линейной комбинации строк матрицы д(Ыz)/дz, подставляя значение д = 1, получим

у дУд(Мгг) = мдУ + Ж = мдУ + Жу = д(МУ) = д± = 0

Таким образом, если производные д^/д^,дКг/дг] вычисляются в предположении линейной однородности У по с^, ^, из уравнений (8) линейной комбинацией с коэффициентами дУ/да, дV/дr| получаем

dVdNi , dVdE (^dVdNi , dV dH\ dp dV dp _

(9)

dci dp ' d] dp I dt'

dV dNi + dV dE _ / ул dV dNi + dci dt dт] dt \i dci dp d] dp J dt dm] dt

I -< dp- dp dm dp I di

Величина в скобках в уравнении (9) имеет размерность Па 1 и будет далее называться коэффициентом сжимаемости

У аут + -УЗЕ (10)

дсг др дг] др

Используя введенное в (10) определение коэффициента сжимаемости, получаем, что комбинация законов сохранения (1) и (5) с коэффициентами дУ/дСг, дУ/д^ дает уравнение

d d V d V d V d V

+ Z d^dlvQ + WdlvJ _ ^ d^ + (11)

Пользуясь законом Дарси (2), получим выражения для Qj, J:

Qi _ - ^2паХг,аК — (Vp - pag) _ -Bj(Vp - pig),

J _ - y,nahaK^(Vp - pag) _ -F(Vp - pg),

(12)

a

где использованы обозначения:

Mi _ K^2naXi,a —, F _ K^2naha —,

a a

_ '}2anaxi,aTf^pa _ '}Zanahai7f^pa

Pi -

Y^an»xi,a ha ijZ

Подставим (12) в (11). При этом получим квазилинейное параболическое уравнение:

др тт^дУ , . дУ

д. - £ ^(в,ур) - -дУ

г

£ ду {*■> - ^(в + - ^(ВД} •

d d V d V

^ ^ ^dlv(®iVp) - dlv(FVp) _

dV („ , _ _ л dV - - (13)

д

Уравнение (13) называется уравнением пьезопроводности. Все коэффициенты уравнения явно либо неявно зависят от давления .

В силу необходимости учёта изменения пористости, коэффициенты линейной комбинации дУ/дСг, дУ/дг] не являются постоянными величинами. Их вычисление производится в блоке фазового равновесия [5].

3.2. Общая схема вычислительного алгоритма Инициализация

Начальные данные задаются в виде полей:

• температуры,

• давления,

• насыщенностей и компонентного состава фаз,

• начальной пористости.

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

Ш!аг по времени

Рассматривается полностью неявная схема для законов баланса количества компонентов и энергии.

• Производится сохранение первичных параметров р, д, 7], чтобы к ним можно было вернуться в случае неудачного шага по времени.

• Итерация начинается с пересчета вторичных величин по значениям первичных р, д, 7]. Для этого в каждой ячейке решается задача фазового равновесия. После определения параметров равновесия вычисляются физические характеристики фаз.

• Значения Ni,E, ( вычисляются с учётом поправок, вносимых источниками скважин-ного типа, химическими реакциями и процессом теплопроводности.

• Решается дискретизированное уравнение пьезопроводности относительно р(т+1) с учётом изменённых параметров Ni,E, (.

• Из уравнений баланса количества компонентов и энергии явным образом находятся ^(т+1), E(m+1) и пересчитываются первичные величины c(m+1), ^(т+1).

• Анализируется сходимость. В случае сходимости итераций производится корректировка величины шага по времени и переход к следующему временному слою.

• Если на каком-либо этапе возникли нефизические значения величин, превышено число итераций либо нарушилось условие Куранта, происходит прерывание итераций с уменьшением шага и повторением процесса с нулевой итерации.

• Переход к началу очередной итерации. 4. Результаты

Для обеспечения адекватности моделирования была проведена предварительная проверка разработанного симулятора и его сравнение с коммерческим симулятором Schlumberger Eclipse 300.

В верификационном расчётном сценарии рассматривался трёхслойный коллектор, центральный слой которого был более проницаемым. Нефтяная фаза состояла из четырёх компонентов, три их которых являлись лёгкими нефтями, а один — тяжелой нефтью. При этом в пластовых условиях в нефти были также растворены два газа. PVT-свойства флюида, закачиваемой воды и газов были взяты из открытых источников. В области имеются две скважины: нагнетающая и откачивающая. Обе описывают приток вещества в соответствии с формулой Писмана. Добывающая скважина осуществляет работу при фиксированном давлении на устье, в то время как на закачивающей указаны как максимальное давление, так и максимальный поток закачиваемой воды. Расчёт вёлся с прогнозом на 10 лет.

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

проницаемых пласта только частично вытеснены. Несмотря на одинаковые характеристики менее проницаемых пластов, фильтрационная картина в них разная за счёт учёта гравитационных сил. Отметим, что в результатах, полученных с помощью Schlumberger Eclipse наблюдается немонотонное поведение решений в областях, где насыщенности близки к критическим значениям.

0.159715 0.69891

Рис. 1. Сравнение картины фильтрационного течения (цветом указана объёмная доля воды). Слева — разработанный симулятор BSD, справа — Schlumberger Eclipse

го <

•е-

CD I tS

ю

QJ

сг

4500 4000 3500 3000 2500 2000 1500 1000 500

Сра внение дебито в нефтч 1

у/

/

■Eclipse BSD

500

1000

1500

2000

2500

3000

3500

4000

Время, дни

Рис. 2. Сравнение дебитов нефти

На рис. 2 приведено сравнение дебитов добытой нефти. Наблюдается практически полное согласование графиков накопленной добычи до момента прорыва центрального слоя, после чего разработанный симулятор BSD демонстрирует меньший темп добычи по сравнению с Eclipse. Разница в накопленной добыче не превосходит 4%. Аналогичные результаты демонстрируются и на графике добычи газа.

Схема модельного расчётного сценария представлена на рис. 3. Расчётная область состоит из двух слоёв, верхний из которых содержит кероген и практически непроницаем. По краям области расположены скважины, закачивающие подогретую воду. В центре располагается добывающая скважина, в которую спущен нагреватель. Вычислительная сетка вблизи нагревателя измельчена. В процессе фильтрации верхний слой нагревается и кероген разлагается, образуя дополнительную пористость и проницаемость (см. рис. 4).

На рис. 5 показана насыщенность воды в момент времени, когда фронт вытеснения только подходит к добывающей скважине. Тем не менее в районе добывающей скважины наблюдается насыщенность, меньшая остаточной водонасыщености. Это связано с тем, что

в этой области происходит испарение воды (и лёгких нефтей) и добыча уже подвижного газа.

Рис. 3. Конфигурация модельной задачи

Рис. 4. Зависимость пористости от времени (в десятках дней) в прискважинных ячейках. Чем позже происходит повышение пористости, тем дальше ячейка находится от нагревателя

2,640е-09

Residual Water Saturation

Рис. 5. Поле насыщенности воды в процессе фильтрации в модельной задаче

5. Заключение

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

Работа выполнена при финансовой поддержке РНФ в рамках гранта №15-11-00015 в лаборатории флюидодинамики и сейсмоакустики МФТИ.

Литература

1. Сивухин Д. В. Общий курс физики. Том 2. Термодинамика и молекулярная физика. М.: Наука, 1979. 552 с.

2. Lake L. W. Enhanced oil recovery. New Jersey: Prentice Hall, 1989.

3. Firoozabadi A. Thermodynamics of hydrocarbon reservoirs. New York: McGraw-Hill, 1999.

4. Болгарский Д. В., Мухачев Г. А., Щукин В. К. Термодинамика и теплопередача. Изд. 2-е, перераб. и доп. М.: Высшая школа, 1975. 495 с.

5. Шевченко А. В.,Цыбулин И. В.,Скалько Ю. И. Оптимизационный алгоритм решения задачи о фазовом равновесии // Труды МФТИ. 2015. Том 7, № 1.

6. Chen Z, Huan G., Ma Y. Computational Methods for Multiphase Flows in Porous Media. Philadelphia: SIAM, 2006.

References

1. Sivukhin D.V. General course of physics. Vol. 2. Thermodynamics and molecular physics. M.: Nauka, 1979. 552 p.

2. Lake L.W. Enhanced oil recovery. New Jersey: Prentice Hall, 1989.

3. Firoozabadi A. Thermodynamics of hydrocarbon reservoirs. New York: McGraw-Hill, 1999.

4. Bolgarskiy D.V., Muhachev G.A., Schukin V.K. Thermodynamics and Heat Transfer. Ed. 2nd, revised. M.: Vischaya Shkola, 1975. 495 p.

5. Shevchenko A.V., Tsybulin I.V., Skalko Y.I. An optimization algorithm for solving phase equilibrium problem. Proceedings of MIPT. 2015. Vol. 7, N 1.

6. Chen Z., Huan G., Ma Y. Computational Methods for Multiphase Flows in Porous Media. Philadelphia: SIAM, 2006.

Поступила в редакцию 26.06.2015.