Software & Systems

Известно, что от общего объема разведанных в России запасов нефти около 11 % приходится на высоковязкие нефти [1]. Основная трудность в их разработке – низкие дебиты скважин. Применение традиционных методов, например заводнения, не позволяет достичь высокого значения коэффициента извлечения нефти ввиду вязкостной неустойчивости фронта вытеснения нефти вытесняющими агентами. Основными методами разработки таких месторождений являются тепловые [2–7], в которых обеспечиваются повышение температуры пласта и снижение вязкости нефти, в том числе и путем внутрипластового горения нефти [3–7].

Один из тепловых методов – с закачкой в пласт воздуха и созданием в нефтенасыщенном пласте подвижного очага низкотемпературного окисления, известный как термогазовое воздейст- вие (ТГВ) [2, 6], считается перспективным для разработки нетрадиционных запасов нефти месторождений баженовской свиты.

Особое внимание в связи с этим необходи- мо уделять вопросам корректного численного моделирования тепловых методов и точности прогноза технологических показателей их примене- ния [6].

Одна из проблем при моделировании – измене- ние смачиваемости пород при нагреве, поскольку при тепловых воздействиях температура плас- та может на 100–500 ºС превосходить началь- ную [8–11]. При этом за счет создания искусственного поля температур смачиваемость изменяется как в пространстве, так и во времени. Все это может кардинально изменить характер фильтрации флюидов и результативность теплового воздействия.

Пока этот вопрос мало изучен. Более того, полученные разными авторами экспериментальные данные весьма противоречивы [8, 11], что определяется рядом объективных факторов [9]. Во-первых, эксперименты проводились по различным методикам, с различными нефтями и в разных условиях. Влияние смачиваемости оценивалось также по-разному: по измерениям краевого угла смачивания, по капиллярному впитыванию, по виду кривых относительных фазовых проницаемостей и т.п. Во-вторых, смачиваемость зависит от многих физических факторов, по-разному реагирующих на повышение температуры [11].

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

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

Упрощенное исследование подобного рода предпринято в работе [10]. Рассматривался вариант, когда за счет нагрева смачиваемость всего пласта скачком менялась на противоположную – с гидрофильной на гидрофобную или наоборот. Показано, что такое резкое изменение смачиваемости существенно меняет результат теплового воздействия при разработке. В ряде случаев это может быть сопоставимо с эффектом уменьшения вязкости нефтей с ростом температуры [10]. Однако равномерный прогрев всего пласта может иметь место лишь в немногих методиках, например, при некоторых способах термолиза in-situ (внутри пласта) горючих сланцев [12], что для России в настоящее время неактуально. По этой причине в данной работе предпринята попытка оценить влияние смачиваемости – переменной и по пространству, и во времени.

Подходы к описанию неизотермической фильтрации

При численном моделировании фильтрации учет переменной смачиваемости, строго говоря, требует учета капиллярных сил, то есть отказа от простой схемы Бакли–Леверетта (БЛ-модель [13]) и перехода к усложненной модели Маскета–Леверетта (МЛ-модель [10]), в изотермическом случае приводящей к уравнению Рапопорта–Лиса [8, 14]. Более сложная модель, названная тепловой моделью Маскета–Леверетта (МЛТ-модель [13]), учитывает явно зависимость капиллярных сил от температуры (коэффициента поверхностного натяжения и краевого угла смачивания).

Соответствующие уравнения в общем виде давно выписаны [8]. Однако для классических коллекторов и изотермического случая капиллярные силы обычно считаются пренебрежимо малыми и рассматривается только схема Бакли–Леверетта. То же предположение делается обычно и при неизотермической фильтрации [8, 13–19]. Если неизотермический характер фильтрации и учитывается, то чаще всего только введением температурной зависимости вязкости [8] и фазовых проницаемостей [8, 17]. Если же капиллярное давление учитывается, то для упрощения считается не зависящим явно от температуры [8, 14], а определяется только величиной насыщенности, как и в уравне- нии Рапопорта–Лиса.

Если учитывать зависимость капиллярных сил от температуры, то в развернутой записи уравнения Рапопорта–Лиса появляются новые члены, содержащие градиент температуры. Это потребует корректировки системы аппроксимации и схемы численного решения, которое может быть получено лишь после экспериментального исследования зависимости от температуры угла смачивания и коэффициента поверхностного натяжения [9, 11].

Полное описание процесса переноса должно включать также учет диффузии, теплопроводности, растворения газов в жидкостях, фазовые переходы, сжимаемость фаз и т.д. [18]. Однако учет этих процессов существенно усложняет задачу гидродинамического моделирования, а главное – затрудняет оценку влияния основных эффектов. Поэтому, следуя большинству работ по многофазной неизотермической фильтрации (обзор таких исследований дается в работе [18]), целесообразно пренебречь диффузионными, капиллярными и неравновесными эффектами (так называемое крупномасштабное приближение). Иными словами, учесть влияние температурного изменения смачиваемости только через изменение фазовых проницаемостей при нагреве. Более того, обычно изменением самого вида кривых пренебрегают; влияние температуры учитывается только через изменение характерных точек на кривых фазовых проницаемостей. Такой подход принят и в данной работе.

Температурная зависимость фазовых проницаемостей

Фазовые проницаемости, в первую очередь, определяются значениями остаточной нефтенасыщенности и водонасыщенности. С ростом температуры коэффициент поверхностного натяжения на границе вода–нефть уменьшается [8, 18]. Во многих экспериментах при температуре 45–50 ºС начинается и уменьшение угла смачивания [8]: темп убывания составляет 3º на 10 ºС. По этим причинам с ростом температуры уменьшается и остаточная нефтенасыщенность (рис. 1) [12].

Аналогичные данные приведены в [8]: остаточная нефтенасыщенность для нефти и пород Арланского месторождения составляет 46, 40, 32, 28 % для температур 10, 24, 50 и 150 °С. В целом темп снижения остаточной нефтенасыщенности по ряду экспериментов составляет примерно 1 % абс. на 10 °С.

Соответственно, по данным значительного числа исследований (см. обзор работ 1970–80-х гг. в [18]), фазовая проницаемость нефти до перехода воды в парообразное состояние с ростом температуры возрастает (рис. 2, 3), а согласно данным работ [8, 20], с увеличением температуры рост относительной фазовой проницаемости (ОФП) для нефти сопровождается одновременным снижением ОФП для воды, то есть увеличением остаточной водонасыщенности, темп этого увеличения составляет примерно 1 % абс. на 11 ºС [8].

В работе [8] при разных вариантах численного моделирования неизотермической фильтрации принят линейный закон изменения остаточной нефтенасыщенности (Sн ост) от температуры (Т):

Sн. ост.(Т) = Sн. ост.(Т0) – 0,0015 (Т – Т0),

Sн. ост.(Т) = 0,28 – 0,0015 (Т – Т0),

Sн. ост.(Т) = 0,20 – 0,0020 (Т – Т0),

где Т0 – начальная температура.

Линейный характер зависимости от температуры [8] принят также и для остаточной водонасыщенности:

Sв. ост. (Т) = 0,2 + 0,001 (Т – Т0).

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

На основании литературных данных были сконструированы модельные кривые ОФП, представ- ленные на рисунке (см. http://www.swsys.ru/uploaded/image/2018-3/2018-3-dop/1.jpg) и зависи- мости от температуры характерных точек на этих кривых (рис. 4).

Математическая модель неизотермической многофазной многокомпонентной фильтрации с химическими реакциями с учетом зависимости ОФП от температуры

Композиционная термохимическая модель фильтрации представляет собой систему уравнений и дополнительных соотношений многокомпонентной трехфазной фильтрации нефти, газа и воды при переменной температуре с неподвижной твердой четвертой фазой и химическими реакциями [21, 22]. Эта модель является одной из самых сложных моделей нефтегазовой подземной гидромеханики [21] и обычно используется для модели- рования различных тепловых методов добычи нефти [3, 7], в том числе внутрипластового горения [22] и термогазового воздействия.

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

В композиционной термохимической модели фильтрации обычно выделяют до четырех фаз [22], при использовании четырех фаз можно выделить от четырех и более компонентов:

-     водная фаза (вода); не смешивается с углеводородными фазами, состоит из одного компонента – воды;

-     жидкая углеводородная фаза (нефть); состоит из смеси углеводородных компонентов, находящихся при данном давлении, температуре и концентрациях других компонентов в жидком состоянии;

-     газообразная углеводородная фаза (газ); состоит из смеси углеводородных компонентов, компонента-кислорода, компонента-воды и других газообразных компонентов, находящихся при данном давлении, температуре и концентрациях других компонентов в жидком состоянии;

-     твердая фаза (кокс); состоит из одного твердого компонента, например кокса или парафина.

Элементарный объем Vb может быть представлен в виде Vb = Vr + Vp, Vp = Vs + Vf, Vf = Vw+ Vo + Vg, где Vr – объем породы (участвует в описании тепловых процессов); Vp – поровый объем; Vf – подвижный объем; Vs – объем твердой фазы; Vα, α = w, o, g – объемы нефтяной, водной и газовой фаз.

Пористость f – доля объема, доступная для заполнения фазами:

.

Подвижная пористость ff – доля объема, доступная для заполнения подвижными фазами:

Насыщенность Sα фазой α (α = w, o, g) – доля объема пористой среды, доступного для подвижных фаз, занимаемая этой фазой:

, α = w, o, g,

Sw + So + Sg = 1.

Насыщенность твердой фазой вводится как , так что,  и приведенные насыщенности флюидами (подвижными фазами)  (α = w, o, g): , α = w, o, g, так что, .

Состояние каждой фазы α (α = w, o, g) задается следующими величинами (неизвестными, опреде- ляемыми в процессе расчета):

-     T = T(t, x, y, z) – температура среды (предполагается, что фазы находятся в термодинамическом равновесии, поэтому температуры всех фаз и породы в точке среды совпадают);

-     Pα = Pα(t, x, y, z) – давление в фазе α (α = w, o, g), в данной работе капиллярным давлением пренебрегается, считается, что давление во всех фазах в элементарном объеме одинаково и равно P = P(t, x, y, z);

-     Sα = Sα (t, x, y, z), α = w, o, g, (t, x, y, z) – насыщенности фазами.

Разделим фазы α, α = w, o, g, s = 1, …, nα, на две группы: подвижные (вода, нефть, газ) α = w, o, g = 1, …, n¢α, n¢α = nα – 1, и твердую α = s = nα. Также разделим компоненты с, с = 1, …, nс, на две группы:

-     с = 1, …, n¢с – компоненты, которые могут быть только в подвижных фазах;

-     с = n¢с + 1, …, nс – компоненты, которые могут быть только в твердой фазе.

Обозначим:

Nc = Nc(t, x, y, z) – молярная плотность компонента c, с = 1, …, n¢с, в подвижном объеме, так что, NcVf = NcfVb – количество компонента c в объе- ме Vb;

Nc = Nc(t, x, y, z) – молярная плотность компонента c, с = n¢с + 1, …, nс, в поровом объеме, так что, NcVp = Nc fVb – количество компонента c в объе- ме Vb.

Распределение компонентов по фазам задается nc ´ na матрицей концентраций xc,a = xc,a(Pa, N), N = (N1, …, ).

Рассмотрим прямолинейно-параллельное течение в горизонтальной плоскости. Капиллярным давлением (как сказано выше) и гравитационными силами будем пренебрегать.

Система уравнений сохранения количества компонентов примет вид

 с = 1, …, n¢с,                                   (1)

, с = n¢с + 1, …, nс,                   (2)

где xa – молярная плотность фазы a;  – отток (в молях) компонента c за счет работы скважин;  – приток (в молях) компонента c за счет химических реакций; wa – скорость фильтрации фазы a.

Скорость фильтрации для фазы a определяется с помощью закона Дарси:

,                                                   (3)

где k – абсолютная проницаемость пласта; kra(Sw, Sg) – относительная проницаемость для фазы a; ma – динамическая вязкость a.

Предполагая выполнение условий локального термодинамического равновесия во всех точках пласта, запишем уравнение сохранения энергии:

 (4)

где Ub – внутренняя энергия элементарного объема; Ha – энтальпия фазы a; Kb – теплопроводность элементарного объема;  – плотность потока энергии за счет химических реакций;  – плотность потока энергии за счет работы скважин,  – плотность оттока энергии за счет потерь тепла через кровлю и подошву.

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

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

Теплопроводность элементарного объема определим как

где Ka  – теплопроводность фазы; Kr – теплопроводность породы.

Стехиометрическое уравнение для химической реакции r записываем как

,

где  – вектор стехиометрических коэффициентов для реагентов Ci в реакции r;  – вектор стехиометрических коэффициентов для продуктов Cj в реакции r.

Скорость реакции r зададим как

где Ar – константа скорости реакции r; Er – энергия активации реакции r;  – веса компонентов, влияющих на реакцию r (для нереагентов nri > 0 в случае катализатора и nri < 0 в случае ингибитора); R – газовая постоянная; , cri – концентрация компонента в фазе, участвующей в реакции; nrp, nrs – показатели степени для учета влияния пористости и доли подвижных фаз на скорость реакции (должны быть неотрицательными), Fr = {1, …, nc}.

Приток компонентов qri,c за счет химических реакций и приток энергии qRe, возникающий в силу химических реакций, определяем как

 

где Hr – энтальпия реакции r; i – номер реагента.

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

Просуммируем уравнения (1) по компонентам с = 1, …, n¢с и перепишем в следующем виде:

,       (5)

где  и .

Поскольку , то .

Так как Cr <<1, то 1 + Cr(P – P0) » 1, следовательно, f » f0.

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

где Dt – размер шага по времени, а верхний индекс n – номер шага по времени.

Запишем дискретный аналог уравнения (5) для прямоугольной блочно-центрированной сетки:

(6)

где Dx, Dy, Dz – размеры ячейки в соответствующих направлениях; Dtn+1 – переменный шаг по времени (определяется в соответствии с рекомендация- ми [21]), нижний индекс i – номер ячейки вдоль оси X; PIa – коэффициент продуктивности скважины по фазе a; Pw – забойное давление. В (6) скорость фильтрации аппроксимируется как

,

где с учетом (3)

, a = w, o, g;

; ;

                          (7)

,  и  взвешиваются вверх по потоку аналогично (7). Для вычисления  используется подход, аналогичный описанному подходу для .

Линеаризуя в уравнении (6) слагаемое  методом Ньютона (пренебрегая при этом производными по температуре и концентрации):

,

получаем неявное уравнение для давления.

Чтобы получить уравнение для определения насыщенности твердой фазой , перепишем уравнение (2) с учетом того, что , , и просуммируем его по c:

, .           (8)

Учитывая правило дифференцирования произведения, получим:

                  (9)

Запишем дискретный аналог уравнения (8) с учетом соотношения (9):

                                                          (10)

Линеаризуя в уравнении (10) слагаемое  методом Ньютона (пренебрегая при этом производными по температуре и концентрации):

,

получаем явное уравнение для насыщенности твердой фазой.

Перепишем уравнение (1) в качестве уравнения сохранения количества подвижных компонентов в подвижных фазах следующим образом:

с = 1, …, n¢с,                                                        (11)

где .

Перепишем уравнение (11) в следующем виде:

,

с = 1, …, n¢с,                                                        (12)

где ,

.

Тогда, согласно правилам дифференцирования,

                 (13)

Запишем дискретный аналог уравнения (12) с учетом соотношения (13):

    (14)

где .

Уравнение (14) является явным уравнением для числа молей компонента c в подвижных фазах Nc,wog.

Мольная доля компонента c в подвижных фазах при этом определяется как

.

Рассмотрим уравнение сохранения энергии (4). Воспользуемся правилом вычисления производной функции многих переменных для кумулятивного члена в левой части уравнения (4) (пренебрегая при этом производными по концентрации):

.                                (15)

Запишем дискретный аналог уравнения (4) с учетом (15):

                                      (16),

где  аналогично для  определяется  как среднее гармоническое между ячейками i и i + 1, потеря тепла через кровлю и подошву  рассчитывается по схеме Винсома и Вестервельда.

Уравнение (16) является неявным уравнением для температуры.

Будем считать c = 1 водным компонентом. Поскольку водный компонент не присутствует в нефтяной фазе и углеводородных компонентов нет в водной фазе, то x1,o º xw,o = 0; xc,w = 0, c Î {2, …, n¢c}; xw,w = 1.

Значение Rg Î [0, 1] (мольная доля газовой фазы в подвижном флюиде) определяется из уравнения

,

где значения констант равновесия компонентов Kc = Kc(P, T) заданы корреляционной формулой

,

коэффициенты которой Ac, Bc, Cc, Dc и Ec являются постоянными для каждого компонента.

Значения Rw (мольная доля воды в подвижном флюиде) и Ro  (мольная доля нефтяной фазы в подвижном флюиде) определяются из соотношений , ,  при заданной зависимости константы равновесия воды от давления и температуры Kw(P, T). Так как  и xw,w = 1, значит, . Для расчета содержания компонента c в нефтяной и газовой фазах xc,o и xc,g воспользуемся соотношениями

,

.

Поскольку , , a = 1, …, n¢a, определим насыщенности жидкими фазами следующим образом:

,  и Sg = 1 –  So – Sw.

Учет зависимости ОФП от температуры при математическом моделировании

При задании ОФП в двухфазных системах задаются приведенные далее функциональные зависимости [21, 22].

ОФП в системе нефть–вода:

,

,

где Swcr – насыщенность связанной водой; Sowcr – насыщенность остаточной нефтью при вытеснении водой; коэффициенты A, B, a, b определяются по результатам интерпретации исследований керна и могут зависеть от пористости или проницаемости.

Аналогично ОФП в системе нефть–газ:

,

,

где Sgcr – насыщенность связанным газом; Sogcr – насыщенность остаточной нефтью при вытесне- нии газом; коэффициенты C, D, g, d определяются по результатам интерпретации исследований керна и могут зависеть от пористости или проницаемости.

Однако часто эффективнее задавать значения ОФП в двухфазных системах в виде двух таблиц от водо- и газонасыщенности.

ОФП нефти при трехфазной фильтрации определяется по второму методу Стоуна:

,

где krocw = krow(Swcr) = krog(Sg = 0) – ОФП по нефти при насыщенности связанной водой в отсутствие газа.

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

Перечисленные в таблице константы можно изменить, задав зависимость от температуры.

Вычисленные масштабированные значения ОФП (со штрихом) должны удовлетворять следующим условиям:

; ;

;

В литературе существуют два варианта масштабирования ОФП – двухточечное и трехточечное. Рассмотрим двухточечное масштабирование ОФП, которое сохраняет значения проницаемостей, заданных в таблицах ОФП, в двух точках (табл. 2) в варианте, изложенном в руководстве [23].

Таблица 1

Константы ОФП

Table 1

Constants of relative phase permeabilities

Величина	Описание величины
Swl	Наименьшее значение Sw в таблице для воды
Swcr	Наибольшее значение Sw в таблице для воды, при котором krw(Sw) = 0
Swu	Наибольшее значение Sw в таблице для воды
Sgl	Наименьшее значение Sg  в таблице для газа
Sgcr	Наибольшее значение Sg в таблице для газа, при котором krg(Sg) = 0
Sgu	Наибольшее значение Sg в таблице для газа
Sowcr	Наибольшее значение So = 1 – Sw – Sgl  в таблице для воды, при котором krow(Sw) = 0
Sogcr	Наибольшее значение So = 1 – Sg – Swl  в таблице для газа, при котором krog(Sg) = 0
krw, max	Максимальное значение функции krw(Sw)
krg, max	Максимальное значение функции krg(Sg)
kro, max	Максимальное значение функции krow(Sw) = 0 и krog(Sg) = 0
krwr	krw(1 – Sowcr – Sgl)
krgr	krg(1 – Sogcr – Swl)
krorg	krog(Sgcr)
krorw	krow(Swcr)

Таблица 2

Сохранение значений табличных ОФП при двухточечном масштабировании

Table 2

Preserving relative phase permeabilities table values during two-point rescaling

ОФП	Точка 1	Точка 2
krw	S¢wcr	S¢wu
krg	S¢gcr	S¢gu
krow	S¢wl	1 – S¢owcr – S¢gl
krog	S¢gl	1 – S¢ogcr  – S¢wl

Для заданных значений немасштабированной насыщенности водой Sw и температуры T масштабированная водонасыщенность в таблице ОФП и значение ОФП для воды вычисляются следующим образом:

.

Для нефти в присутствии воды вычисляются

;

Для заданных значений немасштабированной насыщенности газом Sg и температуры T масштабированная газонасыщенность в таблице ОФП и значение ОФП для газа вычисляются следующим образом:

;

Для нефти в присутствии газа вычисляются

;

Рассмотрим масштабирование ОФП при использовании для расчета ОФП по нефти в случае трехфазного течения второй модели Стоуна (как указано выше). Введем следующие функции:

Тогда ОФП для нефти в трехфазной системе можно записать как

                (17)

. 

Если при расчете ОФП по нефти в трехфазной системе по формуле (17) получается, что kro(Sw, Sg, T) < 0, то полагают kro(Sw, Sg, T) = 0. Так как krow(Swl, T) = krog(Sgl, T) – максимальное значение в таблице ОФП, формула (17) обеспечивает выполнение условий kro(Swl, Sg, T) = krog(Sg, T) и kro(Sw, Sgl, T) = = krow(Sw, T).

Оценка влияния на нефтеотдачу переменной смачиваемости породы пласта-коллектора на примере внутрипластового горения

Для такой оценки была построена одномерная фильтрационная модель со следующими параметрами: число ячеек по осям X´Y´Z – 30´1´1; размеры ячеек по осям X´Y´Z – 6,6666´80´20 м; общие размеры пласта 200´80´20 м; пористость 25 %; горизонтальная проницаемость 700 мД; объемная теплоемкость пласта 2216 кДж/м3/°С; коэффициент теплопроводности 10,13 кДж/м/час/°С. Начальные параметры пласта: пластовая температура 30 °С; пластовое давление 100 атм.; водонасыщенность 20 %; газа нет. При расчете использованы следующие компоненты: О2, СО2, N2, вода и псевдокомпоненты: фракция тяжелой нефти, фракция легкой нефти, кокс. Наличием растворенного в нефти природного газа можно пренебречь, так как высоковязкие нефти обычно содержат небольшое количество растворенного газа. На начальный момент времени углеводородный компонент в пласте представлен только фракцией тяжелой нефти. На рисунке 5 показана зависимость вязкости фрак- ции тяжелой нефти в жидкой фазе от температуры. Вязкость тяжелой нефти несколько завышена, как и в работе [22]. Зависимостью вязкости от дав- ления пренебрегается. Вязкость нефти в начальных пластовых условиях (30 °С и 100 атм.) составляет 1 000 сПз.

Модель химических реакций была взята из работы [22] со следующими реакциями:

-     горение тяжелой нефтяной фракции HEAVY+18,5O2→12CO2+13H2O;

-     горение легкой нефтяной фракции LIGHT+5O2→3CO2+4H2O;

-     пиролиз тяжелой нефтяной фракции HEAVY→2LIGHT+4,67COKE;

-     горение кокса COKE+1,25O2→CO2+0,5H2O,

где HEAVY – фракция тяжелой нефти; LIGHT – фракция легкой нефти; COKE – кокс.

Параметры химических реакций отражены в таблице 3.

Таблица 3

Параметры химических реакций

Table 3

Chemical reaction parameters

Реакция	Скорость реакции, кг-моль/сут	Константа энергии активации, кДж/кг-моль	Энтальпия реакции, кДж/кг-моль
1	0,45Е+06	33300	8,18Е+06
2	0,45Е+06	33300	2,22Е+06
3	0,135Е+06	28800	4,69Е+04
4	0,45Е+06	23400	0,54Е+06

Свойства флюидов взяты из работы [22]. Свойства фракций нефти представлены в таблице 4.

Молекулярная масса кокса принята 13 а.е.м. Плотность кокса при температуре 15 °С и давлении 1 атм. равна 1 200 кг/м3.

Моделировалось сухое внутрипластовое горение. В ячейке с координатами 1´1´1 расположена нагнетательная скважина, а в ячейке с координатами 30´1´1 – добывающая. В нагнетательную скважину закачивается воздух, обогащенный кислородом до 40 % по массе. Температура закачиваемого воздуха на забое нагнетательной скважины – 100 °С. Закачка рабочего агента – 56 тыс. нм3/сут., забойное давление в нагнетательной скважи- не – 150 атм., забойное давление в добывающей скважине – 70 атм. В качестве ограничения на добывающей скважине принята обводненность 98 %, мольная доля кислорода в продукции скважины выше 1 % и газонефтяной фактор более 5 тыс. нм3/м3.

Таблица 4

Свойства фракций нефти

Table 4

Properties of oil fractions

Параметр	Фракция тяжелой нефти	Фракция легкой нефти
Молекулярная масса, а.е.м.	170	44
Критическое давление, атм.	18,5	43,1
Критическая температура, °С	897	735
Плотность, кг/м3	при 1 атм. и 0 °С равна 853	при 70 атм. и 0 °С равна 747
Сжимаемость, 1/атм.	1,0×10-5	2,2×10-4
Коэффициент температурного расширения, 1/К	3,0×10-4	6,2×10-4
Теплота парообразования, кДж/кг	615	213

Рассмотрены два варианта расчетов – без учета зависимости ОФП от температуры и с учетом такой зависимости.

На рисунке 6 приведен полученный в результате моделирования профиль распределения температуры в пласте через 25 лет после начала разработки, на отметке 0 м находится нагнетательная скважина, 200 м – добывающая. Видно, что наличие зависимости ОФП от температуры в расчетах приводит к смещению теплового фронта на 50 метров в сторону добывающей скважины, что составляет 25 % расстояния между скважинами, а следо- вательно, существенно влияет на показатели разработки, в том числе на нефтеотдачу.

На рисунке 7 приведена полученная зависимость накопленной добычи нефти от времени. Добыча прекращена из-за приближения фронта горения к добывающей скважине и прорыва кислорода. Увеличение скорости роста накопленной добычи нефти к концу разработки по обоим вариантам связано с подходом к добывающей скважине так называемого нефтяного вала (области повышенной нефтенасыщенности), расположенного перед фронтом горения. При наличии зависимости ОФП от температуры в процессе разработки добывается 52,79 тыс. т нефти за 34,83 года. При наличии этой зависимости добыча нефти падает до 48,44 тыс. т, а время разработки – до 29,25 года. Таким образом, не учитывая зависимость ОФП от температуры, накопленная добыча нефти завышена на 19,1 % относ., а срок разработки – на 10,0 % относ.

Заключение

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

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

Работа выполнена при поддержке Программы фундаментальных научных исследований государственных академий наук № I.2.П27, заказ № 18-0110_П, НИР № 0065-2018-0110.

Литература

1.     Якуцени В.П., Петрова Ю.Э., Суханов А.А. Нетрадиционные ресурсы углеводородов – резерв для восполнения сырьевой базы нефти и газа в России // Нефтегазовая геология. Теория и практика. 2009. № 4. С. 11–31.

2.     Боксерман А.А., Цуканов А.А., Гришин П.А. Термогазовый метод увеличения нефтеотдачи: состояние и перспективы применения // Эффективные технологии разработки залежей углеводородов: матер. науч.-практич. конф. 2013. Речица, Беларусь. С. 41–44.

3.     Бурже Ж., Сурио П., Комбарну М. Термические методы повышения нефтеотдачи пластов. М.: Недра, 1988. 424 с.

4.     Антониади Д.Г., Гарушев А.Р., Ишханов Б.Г. Настольная книга по термическим методам добычи нефти. Краснодар: Советская Кубань, 2000. 464 с.

5.     Боксерман А.А. Востребованность современных МУН – обязательное условие преодоления падения нефтеотдачи // Нефтяное хозяйство. 2004. № 10. С. 41–46.

6.     Бетелин В.Б., Юдин В.А., Афанаскин И.В., Воль- пин С.Г., Кац Р.М., Королев А.В. Создание отечественного термогидросимулятора – необходимый этап освоения нетрадиционных залежей углеводородов России. М.: Изд-во НИИСИ РАН, 2015. 206 с.

7.     Partha S. Sarathi. In-situ combustion hand book – principles and practices. Final Report. 1998. URL: http://repository.icse.utah. edu/dspace/bitstream/123456789/5336/2/DOE-PC-91008-0374-OSTI_ID-3175-.pdf (дата обращения: 22.04.2018).

8.     Алишаев М.Г., Розенберг М.Д., Теслюк Е.В. Неизотермическая фильтрация при разработке нефтяных месторождений. М.: Недра, 1985. 271 с.

9.     Юдин В.А., Вольпин С.Г. Возможные изменения смачиваемости пород баженовской свиты при термогазовом воздействии // Вестн. кибернетики. 2017. № 2. С. 73–79.

10.  Юдин В.А., Афанаскин И.В. Необходимость учета изменения смачиваемости пород при моделировании тепловых методов добычи нефти // Тр. НИИСИ РАН. 2018. Т. 8. № 1. С. 56–64.

11.   Punase A., Zou A., Elputranto R. How do thermal recovery methods affect wettability alteration? IJPE, 2014, vol. 2014, Art. ID 538021, 9 p.

12.   Brandt A.R. Converting Oil Shale to Liquid Fuels: Energy Inputs and Greenhouse Gas Emissions of the Shell in Situ Conversion Process. Environmental Science & Technology, 2008, vol. 42, no. 19, pp. 7489–7494.

13.  Осокин А.Е. Автомодельные задачи неизотермической двухфазной фильтрации и теплового пограничного слоя: дис. ... к.ф.-м.н. Горно-Алтайск, Новосибирск, 1998. 112 с.

14.   Алексеева К.О., Аксаков А.В. Моделирование задач двухфазной неизотермической фильтрации на основе уравнения Рапопорта–Лиса // Вестн. УГАТУ. 2014. Т. 18. № 1. С. 174–179.

15.   Байгереев Д.Р. Экономичные разностные схемы для за- дачи трехфазной неизотермической фильтрации в переменных глобальное давление, насыщенность, температура: дис. … PhD. Республика Казахстан. Усть-Каменогорск, 2017. 126 с.

16.   Титов А.П. Повышение нефтеотдачи неоднородных по проницаемости пластов в условиях неизотермической фильтрации: дис. … к.т.н. Уфа, 2008. 172 с.

17.  Бочаров О.Б., Осокин А.Е. Численное исследование автомодельных задач неизотермической двухфазной фильтрации // Сибирский журнал индустриальной математики. 2002. Т. 5. № 1. С. 8–19.

18.  Шарафутдинов Р.Ф. Исследование процессов неизотермической фильтрации жидкости и газа с фазовыми переходами: дис. ... к.т.н. Уфа, 1990. 168 с.

19.  Владимиров И.В., Валиев Э.М. Изучение процессов не- изотермической фильтрации в поровых коллекторах при приме- нении тепловых методов в разработке залежей высоковязкой нефти // Проблемы сбора, подготовки и транспорта нефти и нефтепродуктов. 2014. № 2. С. 27–40.

20.  Ермилов О.М., Ремизов В.В., Ширковский А.И., Чугунов Л.С. Физика пласта, добыча и подземное хранение газа. М.: Наука, 1996. 541 с.

21.   Азиз Х., Сеттари Э. Математическое моделирование пластовых систем. М.–Ижевск: Изд-во ИКИ, 2004. 416 с.

22.   Coats K.H. In-situ combustion model. SPE of AIM, 1980, no. 8394, pp. 533–554.

23.   Программа для моделирования процессов разработки нефтегазовых месторождений tNavigator версия 3.0. Технич. руководство. М.: Rock Flow Dynamics, 2009. 782 с.

References

1. Yakutseni V.P., Petrova Yu.E., Sukhanov A.A. Unconventional hydrocarbon resources – a reserve for replenishing the resource base of oil and gas in Russia. Petroleum Geology. Theoretical and Applied Studies. 2009, no. 4, pp. 11–31 (in Russ.).
2. Bokserman A.A., Tsukanov A.A., Grishin P.A. Thermogas method of increasing oil recovery: state and prospects of application. Proc. Sci.-pract. Conf. “Effective Technologies for the Development of Hydrocarbon Deposits”. 2013, Rechitsa, Belarus, pp. 41–44 (in Russ.).
3. Burzhe Zh., Surio P., Kombarnu M. Thermal Methods for Enhanced Oil Recovery. Moscow, Nedra Publ., 1988,
   424 p.
4. Antoniadi D.G., Garushev A.R., Ishkhanov B.G. Handbook on Thermal Methods of Oil Production. Krasnodar, Sovetskaya Kuban Publ., 2000.
5. Bokserman A.A. Demand for Modern EOR is a Prerequisite for Overcoming the Decline in Oil Recovery. Oil Industry. 2004, no. 10, pp. 41–46 (in Russ.).
6. Betelin V.B., Yudin V.A., Afanaskin I.V., Volpin S.G., Kats R.M., Korolev A.V. Creation of a Domestic Thermohydrosimulator is an Indispensible Stage in the Development of Non-Traditional Hydrocarbon Deposits in Russia. Moscow, FGU FNTs NIISI RAN Publ., 2015, 206 p.
7. Partha S. Sarathi. In-situ Combustion Handbook – Principles and Practices. Final Report. BDM Petroleum Technologies. BDM-Oklahoma Inc. Bartlesville. Oklahoma, National Petroleum Technology Office U.S. DEPARTMENT OF ENERGY. Tulsa. Oklahoma. 1998. Available at: http://repository.icse.utah.edu/dspace/bitstream/123456789/5336/2/DOE-PC-91008-0374-OSTI_ID-3175-.pdf/ (accessed April 22, 2018).
8. Alishaev M.G., Rozenberg M.D., Teslyuk E.V. Non-Isothermal Flow in Oil Fields Production. Moscow, Nedra Publ., 1985, 271 p.
9. Yudin V.A., Volpin S.G. Possible changes in the wettability of the rocks of the bazhenov formation during thermogas effects. Proc. in Cybernetics. 2017, no. 2, pp. 73–79 (in Russ.).
10. Yudin V.A., Afanaskin I.V. The need to take into account changes in rock wettability in the modeling of thermal methods of oil production. SRISA RAS Proc. 2018, vol. 8, no. 1, pp. 56–64 (in Russ.).
11. Punase A., Zou A., Elputranto R. How do thermal recovery methods affect wettability alteration? IJPE, 2014,
    vol. 2014, Art. ID 538021, 9 p.
12. Brandt A.R. Converting Oil Shale to Liquid Fuels: Energy Inputs and Greenhouse Gas Emissions of the Shell in Situ Conversion Process. Environmental Science & Technology. 2008, vol. 42, no. 19, pp. 7489–7494.
13. Osokin A.E. Self-Similar Problems of Non-Isothermal Two-Phase Flow and Thermal Boundary Layer. PhD Thesis. Gorno-Altaysk, Novosibirsk, 1998.
14. Alekseeva K.O., Aksakov A.V. Simulation of Two-Phase Non-Isothermal Flow Problems Based on the Rapoport-Lis Equation. Bulletin of USATU. 2014, vol. 18, no. 1, pp. 174–179 (in Russ.).
15. Baygereev D.R. Economical Difference Schemes for the Problem of Three-Phase Non-Isothermal Filtration in Variables Global Pressure, Saturation, Temperature. PhD Thesis. Republic of Kazakhstan. Ust-Kamenogorsk, 2017.
16. Titov A.P. Increase in Oil Recovery of Reservoir Heterogeneous Permeability in Conditions of Nonisothermal Flow. PhD Thesis. Ufa, 2008.
17. Bocharov O.B., Osokin A.E. Numerical Investigation of Self-Similar Problems of Non-Isothermal Two-Phase Flow. J. of Applied and Industrial Mathematics. 2002, vol. 5, no. 1, pp. 8–19 (in Russ.).
18. Sharafutdinov R.F. Investigation of Processes of Liquid and Gas Nonisothermal Flow with Phase Transitions. PhD Thesis. Ufa, 1990.
19. Vladimirov I.V., Valiev E.M. Studying nonisothermal flow processes in porous reservoirs when applying thermal methods in the development of high-viscosity oil deposits. Problems of Collecting, Preparing and Transporting Oil and Oil Products. 2014, no. 2, pp. 27–40 (in Russ.).
20. Ermilov O.M., Remizov V.V., Shirkovskiy A.I., Chugunov L.S. Reservoir Physics, Gas Production and Underground Storage. Moscow, Nauka Publ., 1996, 541 p.
21. Aziz Kh., Settari E. Mathematical Modeling of Reservoir Systems. Moscow–Izhevsk, IKI Publ., 2004, 416 p.
22. Coats K.H. In-situ combustion model. SPE of AIM. 1980, no. 8394, pp. 533–554.
23. The Program for Modeling Processes of Development of Oil and gas Deposits tNavigator 3.0. Technical Manual. Moscow, Rock Flow Dynamics, 2009, 782 p.