CC BY
169
УДК 622.276.654 Дата подачи статьи: 03.05.18
DOI: 10.15827/0236-235X.031.3.482-494 2018. Т. 31. № 3. С. 482-494
Использование метода математического моделирования
для оценки влияния смачиваемости пород на нефтеотдачу на примере внутрипластового горения
И.В. Афанаскин 1, к.т.н.., ведущий научный сотрудник, ivari@ßfanaskin.ru
A.В. Королев 1, к.т.н.., старший научный сотрудник, [email protected]
B.А. Юдин 1, к.ф.-м.н., старший научный сотрудник, [email protected]
1 Федеральный научный центр Научно-исследовательский институт системных исследований РАН, г. Москва, 117218, Россия
В статье рассмотрена значимость одной из проблем корректного численного моделирования тепловых методов добычи нефти - изменение смачиваемости пород при нагреве. При тепловых воздействиях температура пласта может на 100-500 °С превосходить начальную, и за счет создания такого искусственного поля температур смачиваемость пласта изменяется как в пространстве, так и во времени. Это может кардинально изменить характер фильтрации флюидов и результативность теплового воздействия.
Для изучения температурного изменения смачиваемости необходимы достаточно сложные эксперименты, осуществляемые с учетом температурного изменения большого числа влияющих на смачиваемость процессов. До постановки таких сложных работ необходимо оценить значимость этого фактора для прогноза добычи при использовании тепловых методов разработки.
Влияние температурного изменения смачиваемости учитывали через изменение характерных точек на кривых относительных фазовых проницаемостей.
Авторами была рассмотрена задача многофазной многокомпонентной неизотермической фильтрации с химическими реакциями, в которой учитывается зависимость относительных фазовых проницаемостей от температуры. Для решения возникающей системы уравнений предложена оригинальная, неявная по давлению и явная по температуре и концентрациям конечно-разностная схема. Показано, что при моделировании внутрипластового горения пренебрежение изменениями смачиваемости, которые обусловлены созданием в пласте поля температур, переменного по пространству и времени, может приводить к значительным погрешностям в прогнозе накопленной добычи нефти (до 20 % относ.) и времени прорыва кислорода к добывающей скважине (до 10 % относ.). Для учета этого эффекта необходимо экспериментальное изучение изменений смачиваемости (и относительных фазовых проницаемостей) при повышении температуры породы.
Ключевые слова: смачиваемость, зависимость относительных фазовых проницаемостей от температуры, математическое моделирование фильтрации, внутрипластовое горение.
Известно, что от общего объема разведанных в России запасов нефти около 11 % приходится на высоковязкие нефти [1]. Основная трудность в их разработке - низкие дебиты скважин. Применение традиционных методов, например заводнения, не позволяет достичь высокого значения коэффициента извлечения нефти ввиду вязкостной неустойчивости фронта вытеснения нефти вытесняющими агентами. Основными методами разработки таких месторождений являются тепловые [2-7], в которых обеспечиваются повышение температуры пласта и снижение вязкости нефти, в том числе и путем внутрипластового горения нефти [3-7].
Один из тепловых методов - с закачкой в пласт воздуха и созданием в нефтенасыщенном пласте подвижного очага низкотемпературного окисления, известный как термогазовое воздействие (ТГВ) [2, 6], считается перспективным для разработки нетрадиционных запасов нефти месторождений баженовской свиты.
Особое внимание в связи с этим необходимо уделять вопросам корректного численного моделирования тепловых методов и точности прогноза технологических показателей их применения [6].
Одна из проблем при моделировании - изменение смачиваемости пород при нагреве, поскольку при тепловых воздействиях температура пласта может на 100-500 °С превосходить начальную [8-11]. При этом за счет создания искусственного поля температур смачиваемость изменяется как в пространстве, так и во времени. Все это может кардинально изменить характер фильтрации флюидов и результативность теплового воздействия.
Пока этот вопрос мало изучен. Более того, полученные разными авторами экспериментальные данные весьма противоречивы [8, 11], что определяется рядом объективных факторов [9]. Во-первых, эксперименты проводились по различным методикам, с различными нефтями и в разных условиях. Влияние смачиваемости оценивалось также по-разному: по измерениям краевого угла смачивания, по капиллярному впитыванию, по виду кривых относительных фазовых проницаемостей и т.п. Во-вторых, смачиваемость зависит от многих физических факторов, по-разному реагирующих на повышение температуры [11].
Для изучения температурного изменения смачиваемости необходимы достаточно сложные экс-
перименты, осуществляемые с учетом большого числа факторов: начальных и текущих термобарических условий в пласте; состава, плотности и вязкости фильтрующихся углеводородных смесей (нефти, а в случае баженовской свиты - дополнительно и природных битумоидов, продуктов искусственного термолиза, окисления нефти и керогена, газов горения т.п.); наличия многофазного насыщения [9].
Разумеется, до постановки таких сложных работ необходимо оценить значимость этого фактора для результата теплового метода разработки.
Упрощенное исследование подобного рода предпринято в работе [10]. Рассматривался вариант, когда за счет нагрева смачиваемость всего пласта скачком менялась на противоположную - с гидрофильной на гидрофобную или наоборот. Показано, что такое резкое изменение смачиваемости существенно меняет результат теплового воздействия при разработке. В ряде случаев это может быть сопоставимо с эффектом уменьшения вязкости нефтей с ростом температуры [10]. Однако равномерный прогрев всего пласта может иметь место лишь в немногих методиках, например, при некоторых способах термолиза ш-бИи (внутри пласта) горючих сланцев [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].
£ 1.0
Ц'0.0 __Г4-»»
s -7.0 -4.5 -2.0 lof!(^v/a)
Рис. 1. Влияние коэффициента поверхностного натяжения на степень извлечения нефти
Fig. 1. The influence of the surface tension coefficient on oil recovery ratio
а) б) в) г)
Рис. 2. ОФП при разных температурах: для холодной воды и нефти - Kwi, Koi; для горячей воды и нефти - Kw2, K02; для пара и нефти - Kso, Kos. а) Т = 27 и 79,5 °С; б) Т = 23,9 и 204,4 °С; в) и г) Т = 18,3 и 204,4 °С
Fig. 2. Oil relative permeabilities at different temperatures: for cold water and oil - Kwi, Koi; for hot water and oil - Kw2, Ko2; for steam and oil - Kso, Kos а) Т = 27 and 79,5 °С; б) Т = 23,9 and 204,4 °С; в) and г) Т = 18,3 and 204,4 °С
Аналогичные данные приведены в [8]: остаточная нефтенасыщенность для нефти и пород Арлан-ского месторождения составляет 46, 40, 32, 28 % для температур 10, 24, 50 и 150 °С. В целом темп снижения остаточной нефтенасыщенности по ряду экспериментов составляет примерно 1 % абс. на 10 °С.
Соответственно, по данным значительного числа исследований (см. обзор работ 1970-80-х гг. в [1 8]), фазовая проницаемость нефти до перехода воды в парообразное состояние с ростом температуры возрастает (рис. 2, 3), а согласно данным работ [8, 20], с увеличением температуры рост относительной фазовой проницаемости (ОФП) для нефти сопровождается одновременным снижением ОФП для воды, то есть увеличением остаточной во-донасыщенности, темп этого увеличения составляет примерно 1 % абс. на 11 °С [8].
В работе [8] при разных вариантах численного моделирования неизотермической фильтрации
принят линейный закон изменения остаточной нефтенасыщенности (& ост) от температуры (Т): & ост.(Т) = ост.(То) - 0,0015 (Т - Тс), & ост.(Т) = 0,28 - 0,0015 (Т - Т0), & ост.(Т) = 0,20 - 0,0020 (Т - Т0), где Т0 - начальная температура.
Линейный характер зависимости от температуры [8] принят также и для остаточной водонасы-щенности:
5в. ост. (Т) = 0,2 + 0,001 (Т - Т0). После перехода воды в пар фазовые проницаемости испытывают скачок, поскольку пар является несмачивающей фазой, а нефть - смачивающей; после такого фазового перехода фазовая проницаемость нефти снижается, а остаточная нефтенасы-щенность растет [18].
На основании литературных данных были сконструированы модельные кривые ОФП, представленные на рисунке (см. http://www.swsys.ru/up-loaded/image/2018-3/2018-3-dop/1.jpg) и зависимости от температуры характерных точек на этих кривых (рис. 4).
Математическая модель неизотермической многофазной многокомпонентной фильтрации с химическими реакциями с учетом зависимости ОФП от температуры
Композиционная термохимическая модель фильтрации представляет собой систему уравнений и дополнительных соотношений многокомпонентной трехфазной фильтрации нефти, газа и воды при переменной температуре с неподвижной твердой четвертой фазой и химическими реакциями [21, 22]. Эта модель является одной из самых сложных моделей нефтегазовой подземной гидромеханики [21] и обычно используется для моделирования различных тепловых методов добычи нефти [3, 7], в том числе внутрипластового горения [22] и термогазового воздействия.
ст о
1
0.8
^ 0.6
зГ
о
* 0.4 0.2 О
-T=31
m Т-ОП
\ V
\\
\
О 0.2 0.4 0.6 0.8 1 Sw, д.ед.
Рис. 3. ОФП нефти при различных температурах пласта
Fig. 3. Oil relative permeabilities at different reservoir temperatures
100 150
Температура, С
а)
Ь 0,2 -
Остаточная нефтенасыщенность при вытеснении водой Остаточная нефтенасыщенность при вытеснении газом
\
\ V \V
V
—•— -1 1-
0 50 100 150 200 250 300
Температура, С
б)
Рис. 4. Использованные модельные зависимости от температуры характерных точек на кривых
относительных фазовых проницаемостей: остаточных насыщенностей вытесняющего (а) и вытесняемого (б) флюидов
Fig. 4. Used temperature dependencies on the temperature
of critical points on relative phase permeability curves: residual displacing (а) and displaced (б) fluid saturations
Данная модель специально разработана для моделирования таких методов добычи нефти, как закачка пара, закачка горячих жидкостей и газов, внутрипластовое горение, закачка химически активных реагентов.
В композиционной термохимической модели фильтрации обычно выделяют до четырех фаз [22], при использовании четырех фаз можно выделить от четырех и более компонентов:
- водная фаза (вода); не смешивается с углеводородными фазами, состоит из одного компонента - воды;
- жидкая углеводородная фаза (нефть); состоит из смеси углеводородных компонентов, находящихся при данном давлении, температуре и концентрациях других компонентов в жидком состоянии;
- газообразная углеводородная фаза (газ); состоит из смеси углеводородных компонентов, компонента-кислорода, компонента-воды и других газообразных компонентов, находящихся при данном давлении, температуре и концентрациях других компонентов в жидком состоянии;
- твердая фаза (кокс); состоит из одного твердого компонента, например кокса или парафина.
Элементарный объем Уь может быть представлен в виде Уь = Уг + Ур, Ур = У + У/, У/ = У.+ Уо + У&, где Уг - объем породы (участвует в описании тепловых процессов); Ур - поровый объем; У/ - подвижный объем; У5 - объем твердой фазы; Уа, а = V, о, g - объемы нефтяной, водной и газовой фаз.
Пористость ф - доля объема, доступная для заполнения фазами:
ф = ^ = - + -f
- + V
f
V
-
- + V + V
f
Подвижная пористость ф/ - доля объема, доступная для заполнения подвижными фазами:
V к + к - V V ф = _£ = _/-!-1 = ф- ^ =
7 кь кь 4 кь
(
=Ф
1 --
V
vf + V ,
v f s J
V 1 - —-
= Ф
V
P J
Насыщенность 8а фазой а (а = V, о, g) - доля объема пористой среды, доступного для подвижных фаз, занимаемая этой фазой: V
а = V, о, g,
V
s„ = = -
Vf V+V + V
Sw + So + Sg 1.
Насыщенность твердой фазой вводится как фу = ф - Ss j и приведенные
- V
S s = v^ V
так что,
насыщенности флюидами (подвижными фазами)
5а (а = V, о, g): 5а = -Ss ) 5а , а = V, о, g, так что,
+ 5 + 5о + 5 g = 1.
Состояние каждой фазы а (а = V, о, g) задается следующими величинами (неизвестными, определяемыми в процессе расчета):
- Т = Т(/, х, у, 2) - температура среды (предполагается, что фазы находятся в термодинамическом равновесии, поэтому температуры всех фаз и породы в точке среды совпадают);
- Ра = РаЦ, х, у, 2) - давление в фазе а (а = V, о, g), в данной работе капиллярным давлением пре-небрегается, считается, что давление во всех фазах в элементарном объеме одинаково и равно Р = Р^, х, у, 2); _ _
- ^а = & У, х, у, 2), а = V, о, g, = (/, х, у, 2) - насыщенности фазами.
Разделим фазы а, а = V, о, g, 5 = 1, ..., Па, на две группы: подвижные (вода, нефть, газ) а = V, о, g = 1, ., п а, п а = Па - 1, и твердую а = 5 = Па. Также
разделим компоненты с, с = 1,
, на две группы:
- с = 1, ..., п'с - компоненты, которые могут быть только в подвижных фазах;
- с = п'с + 1, ..., пс - компоненты, которые могут быть только в твердой фазе.
0,5
0,4
0,3
J 0,2
0,1
0,0
50
0,4
0,3
0,2
0,1
0,1
0,0
n
Обозначим:
N0 = X, у, 2) - молярная плотность компонента с, с = 1, ..., п'с, в подвижном объеме, так что, NcVf = NCфVb - количество компонента с в объеме ¥ъ,
Nc = х, у, 2) - молярная плотность компонента с, с = п'с + 1, ..., пс, в поровом объеме, так что, NcVp = Nc фVb - количество компонента с в объеме Vb.
Распределение компонентов по фазам задается nc х па матрицей концентраций х^а = хс,а(Ра, №), N = (N1, ..., ^).
Рассмотрим прямолинейно-параллельное течение в горизонтальной плоскости. Капиллярным давлением (как сказано выше) и гравитационными силами будем пренебрегать.
Система уравнений сохранения количества компонентов примет вид
[ф(1 - ^ ) Nc ]:
dx
S ( Xc, a^a Wa)
-qc + qRi,c> с = 1 •••>n'с, d _
-(Ф^ )= q^c , с = n' с + 1, • ..,
dt
(1) (2)
где - молярная плотность фазы а; дс - отток (в молях) компонента c за счет работы скважин; Чы с - приток (в молях) компонента c за счет химических реакций; ^а — скорость фильтрации фазы а.
Скорость фильтрации для фазы а определяется с помощью закона Дарси:
ккга 9Р ,,,
™а =---- , (3)
Иа дх
где к - абсолютная проницаемость пласта; куа^, Sg) - относительная проницаемость для фазы а; Иа - динамическая вязкость а.
Предполагая выполнение условий локального термодинамического равновесия во всех точках пласта, запишем уравнение сохранения энергии:
)=-| X [нл^+КкдТ]+
от дх в дх ^ дх) (4)
+Чке - - Ч1> где иъ - внутренняя энергия элементарного объема; На - энтальпия фазы а; Къ - теплопроводность элементарного объема; дКе - плотность потока
энергии за счет химических реакций; - плотность потока энергии за счет работы скважин, -плотность оттока энергии за счет потерь тепла через кровлю и подошву.
Внутренняя энергия элементарного объема вычисляется из внутренних энергий подвижных фаз, твердой фазы и породы:
u ( p,t ) =
(1 - s. hs[u p,t ) h a (t )-p]| + ss^s ( p,t ) h. (t )
--уьф
+УЬ (1-ф)иг (т).
Энтальпия компонентов в твердой и жидких фазах определяется через их удельную теплоемкость и изменение температуры. Аналогично определяется энтальпия породы. Для определения энтальпии компонентов в газовой фазе добавляется еще теплота парообразования. Энтальпия фазы определяется через концентрацию компонентов в фазе, энтальпию компонентов и их молекулярный вес. Энтальпия воды как функция давления и температуры берется из стандартных таблиц.
Теплопроводность элементарного объема определим как
Кь =ф(1 - ^ + КЛ + К^) +
+фкД +(1 -Ф) к,
где Ка - теплопроводность фазы; К - теплопроводность породы.
Стехиометрическое уравнение для химической реакции г записываем как
Х(5РГС ) ,
где ^ = (^ ) - вектор стехиометрических
коэффициентов для реагентов С в реакции г; БРг = () ^ - вектор стехиометрических коэффициентов для продуктов С) в реакции г. Скорость реакции г зададим как
Л = -ф*-(1 -Б.У - А X
хехр[-Ег/(Л -Т)].П(4
где Аг - константа скорости реакции г; Ег - энергия активации реакции г; (пп.) - веса компонентов, влияющих на реакцию г (для нереагентов п^ > 0 в случае катализатора и п^ < 0 в случае ингибитора); Я - газовая постоянная; си = ф-сп., Cгi -концентрация компонента в фазе, участвующей в реакции; пгр, пга - показатели степени для учета влияния пористости и доли подвижных фаз на скорость реакции (должны быть неотрицательными), ¥г = {1, ..., п}
Приток компонентов д„,с за счет химических реакций и приток энергии дЯе, возникающий в силу химических реакций, определяем как
пг пг
дЛ1,с = Х( 5Рг с - Якгс ) Лг > = X НгЛг ,
г=1 г=1
где Нг - энтальпия реакции г; i - номер реагента.
Для приведенной системы дифференциальных уравнений в частных производных, описывающей многофазную многокомпонентную неизотермическую фильтрацию с химическими реакциями, в
d
n
данной работе предложена оригинальная неявная по давлению и явная по температуре и концентрациям конечно-разностная схема, что позволило использовать достаточно простой вычислительный алгоритм для решения задачи, отличающейся большим количеством неизвестных и физических процессов, в том числе и зависимостью относительных фазовых проницаемостей от температуры.
Просуммируем уравнения (1) по компонентам с = 1, ..., п'с и перепишем в следующем виде: Л Л — П _
|К) = -|Ы- I , (5)
8 О'Л с=1
ÔX
(1 - Ss )(?л )
где N, =(1 -wc = , w + £ w + £ w .
, ^ О О T>w w ^g g
Поскольку ф = ф0 [ 1 + Cr (P - P )] , то
8ф 8P
— = ф0Сг — . 8t 0 r 8t
Так как Cr <<1, то 1 + Cr(P - Po) и 1, следовательно, ф и фо.
Тогда, пользуясь формулой для производной произведения и переходя к дискретному аналогу, получаем
8 1 1
-(фЛ,) и фо o(NT -N^N^r -(рИ+ -pn),
где At - размер шага по времени, а верхний индекс n - номер шага по времени.
Запишем дискретный аналог уравнения (5) для прямоугольной блочно-центрированной сетки:
AxAyAz (фо), {[(N,f -(nJ ] + C (nJ (pn+( - P:)} = = -AyAz [(w,A+ -^wa^]- (6)
- Z {(,»):M[P:+(-(PwГ]}+Z)n,
a=w,o,g c=(
где Ar, Ay, Az - размеры ячейки в соответствующих направлениях; Atn+1 - переменный шаг по времени (определяется в соответствии с рекомендациями [21]), нижний индекс i - номер ячейки вдоль оси X; Pia - коэффициент продуктивности скважины по фазе a; Pw - забойное давление. В (6) скорость фильтрации аппроксимируется как
(wХ0.5 = Z [(,a):+ 0.5 (wa):++0(.5 ] , a=w,o, g
где с учетом (3)
k (k ): P:+( - P: ki+0.5 (kra ),+o.5 P+1 P
( WaL.5 =-"
oo:
n
a )i+0.5
Ax
a = w, o, g;
: 2k,+1k, k,.,, + k
V ^c.a Ji+0.5
(k™ )n+0.5
|( un+i, P+ * p.n
l( kra)n , P+1 < РП -
(7)
(^a )n+0.5 , ( Xc,a )i-
взвешиваются
\n
и (u ,
/,+ 0.5 Vr^i+0.5
вверх по потоку аналогично (7). Для вычисления
\n+1
n+1 a Л+0.5 '
()П_05 используется подход, аналогичный описанному подходу для (^ )
Линеаризуя в уравнении (6) слагаемое (Щ )"+1
методом Ньютона (пренебрегая при этом производными по температуре и концентрации):
п (8Щ
(\n
ÔP
( pn+1 - pn )+( n )n
получаем неявное уравнение для давления.
Чтобы получить уравнение для определения насыщенности твердой фазой , перепишем уравнение (2) с учетом того, что хс = Щ , с е {п/ +1,...,и,}, и просуммируем его по с:
д_ dt
__nc _
(^ss* )= : 4Ri,c, c
= n +1,..., n .
(8)
Учитывая правило дифференцирования произведения, получим:
д_
Ôt
+ Ss ^,Ф0С
Ss )«фс
ÔP
rSs ^ + ^ ^
Ôt s Ôt
(9)
Ôt
Запишем дискретный аналог уравнения (8) с учетом соотношения (9):
(Фо {()П Г )П ]+
)П [(Г-(); ]+с, ()П & )П (+1 - Рп )} =
пс
= I (&,с)П. (10)
с=пС +1
Линеаризуя в уравнении (10) слагаемое (^) П+1
методом Ньютона (пренебрегая при этом производными по температуре и концентрации):
(^s)n+1 ^jp ) n ( Pn1 - P^+(^n
получаем явное уравнение для насыщенности твердой фазой.
Перепишем уравнение (1) в качестве уравнения сохранения количества подвижных компонентов в подвижных фазах следующим образом:
Ô_
Ôt
[Ф fZc (IA +lgSg +lwSw )] =
= Xc,o+ Xc,gÇgwg + le + qRi,c
с = 1, ..., П'с, (11)
где фf =ф(1 - Ss).
n
1
k
i+0.5
c=1
Перепишем уравнение (11) в следующем виде:
[ф fNcwog ] = ~ [wwog ] - + ,
С = 1, ..., n'a,
(12)
где Ncwwog = zc (^S0 + £JSW),
c \ 'Пй o g w y
w = x E w + x E w + x E w .
wog c,^o o c,g vg g c,^w w
Тогда, согласно правилам дифференцирования,
_d_ dt'
-(Ф A, w0g )"ф0 (1 - S» )dNc,wog
+N
(1 - Л ^ ä"-Ф'
dt
dP . dS.
Ht
(13)
Запишем дискретный аналог уравнения (12) с учетом соотношения (13):
^ (i j[i-( ~s. )* ][( NCwg )n -( NCw0g i
+ ( Ncwog )i[l-( Ss )П le, ( Г - P')-( Ncwog h
w ) -( w )
wog /i+0.5 \ wog Л-0.5
( S,)i-( Ц J = -AMz S ^.^[Г - (P f ]} + (i*
ГДе (Wwog )"++0, = S [(XCa)i+o.5 (U ++ 0.5 (Wa):
(14)
n+1 i+ 0.5
Уравнение (14) является явным уравнением для числа молей компонента c в подвижных фазах
Nc,wog.
Мольная доля компонента c в подвижных фазах при этом определяется как
( Zc Г =( Nc, )n+У S ( NCi )n+1.
Рассмотрим уравнение сохранения энергии (4). Воспользуемся правилом вычисления производной функции многих переменных для кумулятивного члена в левой части уравнения (4) (пренебрегая при этом производными по концентрации):
dU _ dU SP dUb дТ
(15)
dt dP dt dT dt Запишем дискретный аналог уравнения (4) с учетом (15):
дР I г' ) ' I дТ
AxAyAz
Atn
dJk]nPn+-Pn)+fêl (Tn+Tn)
-AyAz S [[(U!
-(HaL, (^a)n-0, (Wa)n-+o1.5 ] +
AyAz Ax
[( K )n+0.5 (C1- ry
(16)
-( Kb i T +1- T?)]+(q». )i - S |(Ha^aP/a)n [Pf" -ft )П+1 ]} " ( )П
(„ Г |( Р+1 ^ Р >
где (Н„ I = аналогично для
( а)'+° 1( На)П , Р+1 < Р ,
(Н )П_05 определяется (К )П+05 как среднее гармоническое между ячейками i и /' + 1, потеря тепла через кровлю и подошву ()п рассчитывается по
схеме Винсома и Вестервельда.
Уравнение (16) является неявным уравнением для температуры.
Будем считать c = 1 водным компонентом. Поскольку водный компонент не присутствует в нефтяной фазе и углеводородных компонентов нет в водной фазе, то Х1,о = х„,о = 0; хс^ = 0, c е {2, ...,
п с}; 1.
Значение е [0, 1] (мольная доля газовой фазы в подвижном флюиде) определяется из уравнения
/ ч ^ г (к + к -1)
F (Л ) = Х-^ с + "—I—- = °,
( 8} Xх!-+ ^(Кс + К-1) ,
где значения констант равновесия компонентов Кс = Кс(Р, Т) заданы корреляционной формулой Кс (Р,Т) = (Ас + БС1Р + СсР). ехр[-Бс/(Т-Ес)] ,
коэффициенты которой Ас, Вс, Сс, Ос и Ес являются постоянн^1ми для каждого компонента.
Значения (мольная доля воды в подвижном флюиде) и Яо (мольная доля нефтяной фазы в подвижном флюиде) определяются из соотношений Л = г - Л К , Л = 1 - г - Л (1 - К ) , х = К
V -НТ 8 ^ О М* №
при заданной зависимости константы равновесия воды от давления и температуры К№(Р, Т). Так как
х /х = К (Р,Т) и = 1, значит, х = К .
Для расчета содержания компонента с в нефтяной и газовой фазах хс,о и хс^ воспользуемся соотношениями
С 1 - Zw + Rg ( Kc + Kw - 1)'
zcKc
x =_-_-_
c,g 1 - Zw + Rg (K- + Kw-1).
N R
Поскольку Sa =- "
Ntot = S ( Ncw0g ),
a = 1, ..., n'a, определим насыщенности жидкими фазами следующим образом:
S,„ =-
R„N„,
R Yn
o / y c,
S =-
и Sg — 1 — So — Sw.
Учет зависимости ОФП от температуры при математическом моделировании
При задании ОФП в двухфазных системах задаются приведенные далее функциональные зависимости [21, 22].
X
a=w,o, g
c=1
Z
c
c=1
ОФП в системе нефть-вода:
S - SI
1 Swcr
1 - S„.„ - S...
К (5— ) = А к (5 ) = В
_ <жег —
где Бмсг - насыщенность связанной водой; 8омег -насыщенность остаточной нефтью при вытеснении водой; коэффициенты А, В, а, р определяются по результатам интерпретации исследований керна и могут зависеть от пористости или проницаемости. Аналогично ОФП в системе нефть-газ:
К
( Sg ) = с
( Sg ) = ^
Sg Sgcr
1 - S„„
1 Sogcr Sg
1 — S — S
ogcr gcr
где - насыщенность связанным газом; 8оесг -насыщенность остаточной нефтью при вытеснении газом; коэффициенты С, Б, у, 5 определяются по результатам интерпретации исследований керна и могут зависеть от пористости или проницаемости.
Однако часто эффективнее задавать значения ОФП в двухфазных системах в виде двух таблиц от водо- и газонасыщенности.
ОФП нефти при трехфазной фильтрации определяется по второму методу Стоуна:
+ к„ ()
ko, ( Sw, Sg ) = krc
X S. )
+ k.
X Sg )
- kw ( Sw )-krg ( Sg )
где кГоск = кго^сг) = кГог(8г = 0) - ОФП по нефти при насыщенности связанной водой в отсутствие газа.
При задании ОФП так или иначе (с помощью таблиц или функциональных зависимостей) должны быть указаны константы, приведенные в таблице 1.
Перечисленные в таблице константы можно изменить, задав зависимость от температуры.
Вычисленные масштабированные значения ОФП (со штрихом) должны удовлетворять следующим условиям:
у ^ с/ . 5/ <! _ 5
Su *1 - SWi
+ s'< 1 ; S,
+ Sgcr + S'w, < 1.
В литературе существуют два варианта масштабирования ОФП - двухточечное и трехточечное. Рассмотрим двухточечное масштабирование ОФП, которое сохраняет значения проницаемостей, заданных в таблицах ОФП, в двух точках (табл. 2) в варианте, изложенном в руководстве [23].
Константы ОФП Constants of relative phase permeabilities
Таблица 1 Table 1
Величина Описание величины
Swl Наименьшее значение Sw в таблице для воды
Swcr Наибольшее значение Sw в таблице для воды, при котором krWSw) = 0
Swu Наибольшее значение Sw в таблице для воды
Sgl Наименьшее значение Sg в таблице для газа
Sgcr Наибольшее значение Sg в таблице для газа, при котором krg(Sg) = 0
Sgu Наибольшее значение Sg в таблице для газа
Sowcr Наибольшее значение So = 1 - Sw - Sgl в таблице для воды, при котором krow(Sw ) = 0
Sogcr Наибольшее значение So = 1 - Sg - Sw в таблице для газа, при котором 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 g
krow S wl 1 — S owcr — S gl
krog S'gl 1 — S ogcr — S wl
Для заданных значений немасштабированной насыщенности водой Sw и температуры Т масштабированная водонасыщенность в таблице ОФП и значение ОФП для воды вычисляются следующим образом:
, 5—ег (Т)> 5т (Т); 5—ег , 5— < 5—ег (Т); (,Т) = 1 5^, 5— > 5„ (Т);
^ ^ _ ^ег (Т)] "
к' (Т)
к!/ (5— ,Т)= к^-1 к„ (5((5— ,Т)) .
У
s
rw.iiiax
Для нефти в присутствии воды вычисляются
si, ( Sw ,T ) =
Sw, S'wl(T )> 1 - S/wr (T )-SgjT ) ; S„,, Sw < S'wl (T ) ;
1-Swr - Sgl, Sw > 1-SOwr (T )-S!gl(T ) ;
1 + [Sw Swi (T)] 1-S^r (TyS'j gT)-s;, (T) С ( S» ,T )=^p^ krow (s/( Sw ,T )) ;
Slow (Sw,T) =
Sw, S'w, (T)> Sm (T); S„, Sw < S'w, (T ) ; Swu, Sw > S'wu (T) ;
S№l +[ Sw - sw, (T )]■
Swu Swl
х (т)- ^ (т)
Для заданных значений немасштабированной насыщенности газом и температуры Т масштабированная газонасыщенность в таблице ОФП и значение ОФП для газа вычисляются следующим образом:
к/ (, )= к^ (5/ З^,))к
Sg (Sg,T) =
Sg, С (T)> Sgu T;
Sgcr, Sg < С (T ) ;
Sgu, Sg > Su (t ) ; Sgcr +[Sg - Sgcr (T)]-
Sgu Sgcr
'sgu (t)-Sg„ (t) Для нефти в присутствии газа вычисляются
К> Sgi(T)* 1-С(T)-Sl,(T);
S' (S ,T) =
°g V g f
Sgi, Sg<S'gi(T);
1 - S - S ,, S > 1-S' (T )-S', (T );
Ogcr wl> g Ogcr \ ! wl \ / '
Sgl + [Sg - S'gl(T)]-
SOgcr (T )-S'gl(T )-Sw,(T )
Kg ( S, )=kog S ( S, ));
Sg , к, (t )> Sg» (t ) Sgl, Sg < Sgl (t) Sg» , Sg > Sg» (t )
Sig ( Sg,t )=
Sg +[ Sg - Sg (t )]:
s„„ - s„,
gl 'Lsg /Jsg„(t)-sgl(t)
Рассмотрим масштабирование ОФП при использовании для расчета ОФП по нефти в случае трехфазного течения второй модели Стоуна (как указано выше). Введем следующие функции:
c(T ) =
= (T )
.(T )'
Cmax (T )'
k/-g ,max (T }
если заданы кге (Т ) и кп (Т ); если задано только к^ (Т); если задано только к (Т);
rw,max
1 иначе;
ß(T ) =
(Т )
krg ,max (Т )
(T)
, если заданы krg (T) и kw (T);
, если задано только kTO (T) ;
k_
(T )'
если задано только
Ks (т) ;
1 иначе.
Тогда ОФП для нефти в трехфазной системе можно записать как
k/o ( Sw, Sg, T )= k/ow ( Swl (T ) ,T )
Kw ( Sw ,T )
X
Kw ( sw, (T ) ,T )
Kog ( Sg, T ) + Kog ( Sgl (T ) , T )
+ a(T ) k'w ( Sw T )
ß(T ) Kg ( Sg,T )
(1 S)
-[a(T)^ (^,T) + ß(T)Kg (Sg,T)]
Если при расчете ОФП по нефти в трехфазной системе по формуле (17) получается, что km(Sw, Sg, T) < 0, то полагают km(Sw, Sg, T) = 0. Так как kmw(Swi, T) = krog(Sgi, T) - максимальное значение в таблице ОФП, формула (17) обеспечивает выполнение условий kro(Swi, Sg, T) = krog(Sg, T) и kro(Sw, Sgi, T) =
krow(Sw, T)•
Оценка влияния на нефтеотдачу переменной смачиваемости породы пласта-коллектора на примере внутрипластового горения
Для такой оценки была построена одномерная фильтрационная модель со следующими параметрами: число ячеек по осям XxYxZ - 30x1x1; размеры ячеек по осям XxYxZ - 6,6666x80x20 м; общие размеры пласта 200x80x20 м; пористость 25 %; горизонтальная проницаемость 700 мД; объемная теплоемкость пласта 2216 кДж/м3/°С; коэффициент теплопроводности 10,13 кДж/м/час/°С. Начальные параметры пласта: пластовая температура 30 °С; пластовое давление 100 атм.; водонасы-щенность 20 %; газа нет. При расчете использованы следующие компоненты: О2, СО2, N2, вода и псевдокомпоненты: фракция тяжелой нефти, фракция легкой нефти, кокс. Наличием растворенного в нефти природного газа можно пренебречь, так как высоковязкие нефти обычно содержат небольшое количество растворенного газа. На начальный момент времени углеводородный компонент в пласте представлен только фракцией тяжелой нефти. На рисунке 5 показана зависимость вязкости фракции тяжелой нефти в жидкой фазе от температуры. Вязкость тяжелой нефти несколько завышена, как и в работе [22]. Зависимостью вязкости от дав-
k
rg, max
ro.max
X
rg
1-Sgcr - Sgl- S
k
rw.max
„1000
r
900
n
h 800
0
ГЛ 700
IK
600
500
400
300
200
100
0
• •—»
0 25 50 75
100 125 150 175 200 225 250 275 300 Температура, 0С
Рис. 5. Зависимость вязкости фракции тяжелой нефти в жидкой фазе от температуры
Fig. 5. Temperature dependence of heavy fraction viscosity in a liquid phase
ления пренебрегается. Вязкость нефти в начальных пластовых условиях (30 °С и 100 атм.) составляет 1 000 сПз.
Модель химических реакций была взята из работы [22] со следующими реакциями:
- горение тяжелой нефтяной фракции HEAVY+18,5O2^12CO2+13H2O;
- горение легкой нефтяной фракции LIGHT+502^3 CO2+4H2O;
- пиролиз тяжелой нефтяной фракции HEAVY^2LIGHT+4,67 COKE;
- горение кокса 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.
Моделировалось сухое внутрипластовое горение. В ячейке с координатами 1x1x1 расположена нагнетательная скважина, а в ячейке с координатами 30x1x1 - добывающая. В нагнетательную скважину закачивается воздух, обогащенный кислородом до 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,010-5 2,2-10-4
Коэффициент температурного расширения, 1/К 3,010-4 6,2-10-4
Теплота парообразования, кДж/кг 615 213
Рассмотрены два варианта расчетов - без учета зависимости ОФП от температуры и с учетом такой зависимости.
На рисунке 6 приведен полученный в результате моделирования профиль распределения температуры в пласте через 25 лет после начала разработки, на отметке 0 м находится нагнетательная скважина, 200 м - добывающая. Видно, что наличие зависимости ОФП от температуры в расчетах приводит к смещению теплового фронта на 50 метров в сторону добывающей скважины, что составляет 25 % расстояния между скважинами, а следо-
О 250
¡5 150
*
Sr» V
1 \
г т J
150 175 200 Расстояние, м
• — Без учета зависимости ОФП от температуры — С учетом зависимости ОФП от температуры
Рис. 6. Профиль распределения температуры в пласте через 25 лет после начала разработки
Fig. 6. Reservoir temperature profile after 25 years ofproduction
200
100
50
0
0
25
50
75
100 125
вательно, существенно влияет на показатели разработки, в том числе на нефтеотдачу.
На рисунке 7 приведена полученная зависимость накопленной добычи нефти от времени. Добыча прекращена из-за приближения фронта горения к добывающей скважине и прорыва кислорода. Увеличение скорости роста накопленной добычи нефти к концу разработки по обоим вариантам связано с подходом к добывающей скважине так называемого нефтяного вала (области повышенной нефтенасыщенности), расположенного перед фронтом горения. При наличии зависимости ОФП от температуры в процессе разработки добывается 52,79 тыс. т нефти за 34,83 года. При наличии этой зависимости добыча нефти падает до 48,44 тыс. т, а время разработки - до 29,25 года. Таким образом, не учитывая зависимость ОФП от температуры, накопленная добыча нефти завышена на 19,1 % относ., а срок разработки - на 10,0 % относ.
Заключение
Для системы дифференциальных уравнений в частных производных, описывающей многофазную многокомпонентную неизотермическую фильтрацию с химическими реакциями, предложена оригинальная, неявная по давлению и явная по температуре и концентрациям конечно-разностная схема. Это позволяет использовать достаточно простой вычислительный алгоритм для решения задачи с большим количеством неизвестных и физических процессов. Кроме обычно учитываемых при решении такого рода задач явлений и процессов, рассматриваемая модель также учитывает зависимость относительных фазовых проницаемо-стей от температуры. Для обеспечения устойчивости численной схемы необходимо следить, чтобы изменения неизвестных (давление, температура, концентрации) за временной шаг не превышали заранее заданных предельных значений.
При моделировании внутрипластового горения пренебрежение изменениями смачиваемости, обусловленными созданием в пласте поля температур, переменного по пространству и времени, может приводить к значительным погрешностям в прогнозе накопленной добычи нефти (до 20 % относ.) и времени разработки (до 10 % относ.). Учет этого эффекта может быть осуществлен только на основании результатов экспериментального изучения изменений смачиваемости (и относительной фазовой проницаемости) при повышении температуры.
Работа выполнена при поддержке Программы фундаментальных научных исследований государственных академий наук № 1.2.П27, заказ № 18-0110 П, НИР № 0065-2018-0110.
Литература
1. Якуцени В.П., Петрова Ю.Э., Суханов А.А. Нетрадиционные ресурсы углеводородов - резерв для восполнения сырьевой базы нефти и газа в России // Нефтегазовая геология. Теория и практика. 2009. № 4. С. 11-31.
2. Боксерман А.А., Цуканов А.А., Гришин П.А. Термогазовый метод увеличения нефтеотдачи: состояние и перспективы применения // Эффективные технологии разработки залежей углеводородов: матер. науч.-практич. конф. 2013. Речица, Беларусь. С. 41-44.
3. Бурже Ж., Сурио П., Комбарну М. Термические методы повышения нефтеотдачи пластов. М.: Недра, 1988. 424 с.
4. Антониади Д.Г., Гарушев А.Р., Ишханов Б.Г. Настольная книга по термическим методам добычи нефти. Краснодар: Советская Кубань, 2000. 464 с.
5. Боксерман А.А. Востребованность современных МУН - обязательное условие преодоления падения нефтеотдачи // Нефтяное хозяйство. 2004. N° 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/D0E-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. Байгереев Д.Р. Экономичные разностные схемы для за-
Рис. 7. Зависимость накопленной добычи нефти от времени
Fig. 7. Time dependence of accumulated oil production
дачи трехфазной неизотермической фильтрации в переменных глобальное давление, насыщенность, температура: дис. ... 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 с.
Software & Systems Received 03.05.18
DOI: 10.15827/0236-235X.031.3.482-494 2018, vol. 31, no. 3, pp. 482-494
A mathematical model approach for evaluation of oil recovery changes caused by wettability
during in-situ combustion
I.V. Afanaskin l, Ph.D. (Engineering), Leading Researcher, [email protected] A. V. Korolev l, Ph.D. (Engineering), Senior Researcher, [email protected] V.A. Yudin l, Ph.D. (Physics and Mathematics), Senior Researcher, [email protected]
1 Federal State Institution "Scientific Research Institute for System Analysis of the Russian Academy of Sciences" (SRISA RAS), Moscow, 117218, Russian Federation
Abstract. The paper considers a problem arising in simulation of EOR thermal methods that is wettability changes influenced by formation heating. EOR thermal methods may make a formation temperature rise up to 100-500 °C above the initial value. Wettability may change in time and space as a result of such artificial temperature distribution. Thus, the fluid flow in formation and the result of applying EOR thermal methods may change dramatically.
Investigation of wettability temperature changes requires complicated experiments that take into account temperature changes of a big number of physical processes effecting wettability. Before conducting such experiments, it is mandatory to estimate the importance of such factor in forecasting EOR thermal methods result using computer simulation.
The article took into account wettability temperature changes by implementing temperature dependence of critical points on phase relative permeability curves, both for oil and water.
The authors consider the problem of multi-phase multi-component non-isothermal filtration through with chemical reactions that took into account a temperature dependence of relative phase permeabilities. They proposed an original finite-difference scheme that is implicit on pressure and explicit on temperature and concentrations. The calculations showed that if wet-tability changes due to temperature profile variability in formation space and time are neglected, considerable errors may occur in estimates of cumulative production and time of oxygen breakthrough in predicting in-situ combustion results by computer simulation. These errors may be as much as 20 % in cumulative production and 10 % in breakthrough time. To take into account this effect, it is necessary to study the changes in wettability (and relative permeability) experimentally with increasing formation temperature.
Keywords: wettability, temperature dependence of relative phase permeabilities, simulation of flow through porous media, in-situ combustion.
Acknowledgements. The article has been supported by the Program of Fundamental Scientific Research of State Academies of Sciences no. I.2.P27, order no. 18-0110 P, Research no. 0065-2018-0110.
