CC BY
20
УДК 519.63
В. Е. Родионов, И. В. Цыбулин, А. А. Карпаев
Московский физико-технический институт (государственный университет)
Численное моделирование волны внутрипластового горения при двухфазном фильтрационном течении
Рассматривается задача распространения волны внутрипластового горения в рамках модели неизотермической двухфазной фильтрации. Предложены критерии, позволяющие предсказать возможность реализации различных режимов распространения, и алгоритм расчета скорости волны. Проведено численное моделирование реагирующей двухфазной фильтрации и сравнение результатов расчетов с теоретическими оценками.
Ключевые слова: неизотермическая двухфазная фильтрация, внутрипластовое горение, численное моделирование.
V. E. Rodionov, I. V. Tsybulin, A. A. Karpaev
Moscow Institute of Physics and Technology (State University)
Numerical simulation of an in-situ combustion wave for two-phase flow through porous media
The problem of an in-situ combustion wave propagation is studied. The model of a non-isothermal two-phase flow through porous media is considered. The criteria allowing us to predict the feasibility of different propagation regimes and then algorithm for calculating the wave speed are proposed. Numerical simulation of reacting two-phase flow is performed. Comparison of the calculation results with the theoretical estimates is implemented.
Key words: non-isothermal two-phase flow, in-situ combustion, numerical simulation.
1. Введение
В настоящее время повышенный интерес проявляется к технологиям разработки месторождений Баженовской свиты. Этот интерес обусловлен тем, что предполагаемые запасы углеводородов в этих месторождениях оцениваются как 1 трлн тонн. Одной из таких перспективных технологий является термогазовый метод, основанный на закачке в пласт кислородосодержащей смеси и ее трансформации в высокоэффективный смешивающийся с пластовой нефтью вытесняющий агент за счет самопроизвольных внутрипластовых окислительных процессов. Применительно к методам повышения нефтеотдачи в основу метода внутрипластового горения положен процесс сжигания части нефти или кокса/керогена, содержащейся в пористой среде, приводящий, в частности, к нагреву и увеличению подвижности несгоревших фракций. Образовавшиеся при этом продукты реакции вытесняют нефть из пласта. Эффективности метода способствуют довольно высокие начальные пластовые температуры, лежащие в диапазоне 80-140 С. Предполагается, что прогрев нед-ренируемых зон до температур 250-350 С позволит извлечь из них значительную долю нефти [1]. В ряде случаев этот метод обладает большим потенциалом по сравнению с другими способами повышения нефтеотдачи.
В данной работе рассматриваются некоторые вопросы распространения волн внутри-пластового горения. Мы не будем рассматривать сопутствующие горению процессы испарения/конденсации воды и легких фракций нефти, крекинга и пиролиза керогена и др. [2], а подробно сосредоточимся на основном процессе - распространении волны горения. Цель
работы - провести аналитическое исследование и численное моделирование внутрипласто-вого горения при двухфазном течении в пористом слое. Случай однофазного фильтрационного течения газообразного окислителя (или смеси, содержащей окислитель) подробно исследовался многими авторами и по разным поводам. Мы ограничимся ссылкой на монографию [3]. Применительно к проблеме внутрипластового горения аналогичные задачи рассматривались в работах [4] - [8].
Работа организована следующим образом: в п. 2 формулируется постановка задачи, в п. 3 рассматривается автомодельное решение типа «бегущая волна» и проводятся качественные оценки скорости волны внутрипластового горения, в пп. 4 и 5 исследуются два режима распространения волны горения, п. 6 посвящен численному моделированию внут-рипластового горения и сравнению его результатов с теоретическими оценками.
2. Постановка задачи
Рассматривается пористый слой, заполненный частично жидкой нефтью (или неподвижным коксом/керогеном), частично газовой фазой. Через левую границу в слой нагнетается газ, содержащий окислитель, который при достаточно высоких температурах может вступать в реакцию с нефтью. На левой границе слоя тем или иным способом инициируется реакция горения, которая затем распространяется слева направо. В зависимости от параметров задачи (теплофизических свойств фаз и скелета, начальной насыщенности слоя нефтью, концентрации окислителя в нагнетаемом газе и др.) могут реализоваться различные сценарии распространения волн: вытеснения, горения и тепловой. Мы во всех случаях полагаем, что волна вытеснения является наиболее быстрой и режим распространения внутрипластового горения определяется соотношением между скоростями волны горения и тепловой волны. Режим, в котором тепловая волна распространяется быстрее волны горения, называется Reaction-trailing front (RTF), в противном случае режим называется Reaction-leading front (RLF). В первом случае горячие продукты реакции опережают зону горения и фронт реакции распространяется по предварительно прогретой среде. Во втором случае волна горения распространяется по холодному начальному фону. Качественное распределение температуры, насыщенности жидкой фазы и концентрации окислителя для обоих режимов показаны на рис. 1. Рассмотрим следующую постановку задачи. Известно исходное состояние пористого слоя, через левую границу которого с некоторым темпом закачивается флюид, содержащий окислитель. В слое инициируется экзотермическая реакция окисления нефти, которая распространяется по слою, формируя волну горения. Следует понять, в каком режиме и с какой скоростью будет распространяться эта волна, какая температура будет достигаться в зоне реакции, какая часть нефти или окислителя (в зависимости от режима распространения волны) вступит в реакцию.
Указанную задачу рассмотрим в одномерной постановке в рамках модели двухфазной фильтрации. Жидкая фаза содержит один компонент - нефть, а газовая фаза содержит окислитель (кислород) и инертные газы (азот, продукты реакции: CO, CO2, пары воды и др.). Примем, что парциальные объемы компонентов газовой фазы не зависят от давления и температуры (точнее этой зависимостью можно пренебречь), теплоемкости жидкой и газовой фаз постоянны и равны между собой. В этих предположениях мы объединяем инертные газы в один компонент газовой фазы. Окислитель представляет второй компонент. Примем, что в процессе реакции сохраняется объем, т.е. объем продуктов реакции равен объему веществ, вступающих в реакцию. Эти предположения не являются принципиальными, но заметно упрощают анализ задачи. Отметим, что изменение объема при реакции (в том числе из-за изменения температуры) может при определенных условиях привести к развитию неустойчивости и автоколебательному режиму распространения волны горения даже в одномерном случае [9]. Мы пренебрегаем теплопотерями через подошву и кровлю слоя, а также диффузией компонентов.
Волна горения 0>0 Т=То Тепловая волна
I
Т=ъ )
У=У1 8=8о
у
Э=81 I у=0
Тепловая волна^^ Т=Т1 Волна горения 0>0 Т=Т0
1
У=У1 г 3=8о
Б=0 \ У=Уо
Рис. 1. Распределение температуры (синяя линия), насыщенности жидкой фазы (черная линия) и концентрации окислителя (красная линия) при распространении волны горения в ИХЕ-режиме (слева) и КЬР-режиме (справа)
В указанных предположениях одномерное фильтрационное течение описывается системой уравнений:
т^г + _ -а\Л,
дЬ дх
тд+ дудЬ _ ^^ (1)
д д х
тд(1 - У)(1 - 8) + д(1 - Уд _ _
т-^--1--^- _ а3Л
д д х
Здесь 5 - насыщенность жидкой фазой, у - объемная концентрация окислителя в газовой фазе, Шь- скорости фильтрации жидкой и газовой фаз, т - пористость, а\,а2,а3 -объемные стехиометрические коэффициенты.
Углеводороды обладают способностью вступать в экзотермические реакции с кислородом, что может быть использовано для получения тепла непосредственно в нефтяном пласте. Внутрипластовое горение нефти представляет собой гетерогенный процесс, в котором реагирующие вещества находятся в разных агрегатных состояниях, а взаимодействие может происходить только на поверхности химически активных фаз. При окислении углеводородов можно получить широкий спектр продуктов реакции. В некоторых случаях углеводородные цепи разрушаются полностью, углерод связывается с образованием СО2 и/или СО, а водород - с образованием Н2О, в этом случае говорят о сгорании углеводорода. В других случаях углеводородная цепочка сохраняется, но в молекулу внедряются атомы кислорода [10]. Скорость химической реакции примем равной
л _Лм) е -1,
с /( 5, у) „ ,
где с - температура активации,--предэкспоненциальный множитель, ¿(в, у) описывает зависимость скорости реакции от объемных долей фаз в поровом пространстве и концентрации окислителя.
Условие сохранения объема при реакции означает, что а3 _ а\ + а2. Суммируя урав-дШ
нения (1), получаем —— _ 0, где Ш _ + Шс - полная скорость фильтрации. Мы
д х
не учитываем капиллярный скачок давления между фазами и в соответствии с моделью Баклея-Леверетта, считаем, что
_ рШ, _ (1 -
где (р = p(s) - функция Баклея (доля жидкой фазы в потоке), определяемая относительными проницаемостями и вязкостями фаз.
Уравнение баланса энергии в сделанных предположениях запишется в виде
d ( \ dv д ( дТ\
— ¡mcT + mqs + (1 - m)c8T) - + — i WcT + qWL - \—j =0. (2)
Здесь T - температура (одинаковая в обеих фазах и скелете), р - давление, с - объемные теплоемкости жидкой и газовой фаз (считаются равными), cs - объемная теплоемкость скелета, q - объемная калорийность нефти, Л = const - эффективный (усредненный по фазам и скелету) коэффициент теплопроводности.
Перед волной горения (после прохождения волны вытеснения) насыщенность жидкой фазы составляет s0. На левой границе при х = 0 задается темп закачки (скорость фильтрации) газа Wg и концентрация окислителя в нем у\. В случае RTF-режима распространения волны горения задается температура закачиваемого газа Т\, в случае RLF-режима задается исходная температура пористого слоя То. Считается, что при высоких температурах нефть и окислитель реагируют достаточно быстро, так что в RTF-режиме перед фронтом реакции отсутствует окислитель, а в RLF-режиме за фронтом реакции отсутствует нефть.
3. Автомодельное решение «бегущая волна»
Для анализа режимов распространения волны горения вдали от границ (когда их влиянием можно пренебречь) рассмотрим решения уравнений (1) - (2), имеющие вид бегущих волн: Т = Т(z), s = s(z), у = y(z), где z = х — Dt — автомодельная переменная, D > 0 — постоянная скорость бегущей волны. В результате мы переходим в движущуюся вместе с фронтом систему координат, в которой температура, насыщенность и концентрация не меняются со временем. В соответствии с классической теорией горения [11] величина скорости волны D не произвольна, а является решением некоторой нелинейной задачи на собственные значения. Для существования решения вида «бегущая волна» необходимо, чтобы при низких температурах скорость реакции обращалась в ноль. Для этого вводится температура зажигания реакции Т*, ниже которой реакция не протекает.
Интегрирование уравнений баланса объема дает соотношения
W (z) = const,
—m D(a2s — a1(1 — s)y) + — ai(1 — ф)у) = const (3)
и уравнение для структуры фронта волны горения
— mDs + ) = —aiR. (4)
Интегрирование уравнения энергии дает
dT
—D(mcT + msq + (1 — m)csT) + Dp + W (cT + qip) — A— = const. (5)
В дальнейшем считаем, что перепад давления на фронте волны горения мал и пренебрегаем слагаемым Dp в (5). Введем безразмерную скорость фильтрации k = W/(mD), стехиомет-рическое отношение объемов вступающих в реакцию окислителя и горючего А = a2/(ii и
(1 — Cg
коэффициент к* = 1 +--- > 1. Величина к* имеет смысл отношения скорости филь-
тс
трации к скорости тепловой волны [12]. В этих обозначениях уравнения (3) - (5) примут вид
As — (1 — s)y — к(А<р — (1 — ф)у) = const, (6)
ds a1R/mD (7)
dz = 1 — kf' , ( )
Л йТ
(к * — к)сТ - а(к р - 8) +--= еоп81. (8)
ти аг
Естественно считать, что величины Т, в, у постоянны вдали от фронта:
Т (-то)=Ть «(—то) = «1, у(—ю) = уг, Т (+то) = То, в(+то) = во, у(+то) = уо.
В дальнейшем индексом «0» будем помечать (предельные) величины перед фронтом волны горения (справа), индексом «1» — за фронтом (слева).
Из (6), (8) получаем соотношения, связывающие величины перед и за волной горения:
Авг - (1 - «г)уг - к(А<рг - (1 - рг)ш) = ^о - (1 - «о)Уо - к(Аро - (1 - Ро)Уо), (9)
( К - к)с(Тг - То) - дк(рг - ро) + Фг - во) = 0. (10)
Как уже отмечалось, в рамках классической теории горения скорость волны И не произвольна, а является решением некоторой нелинейной задачи на собственные значения. Применительно к рассматриваемой здесь проблеме скорость волны И должна быть выбрана так, чтобы на фазовой плоскости ( в, Т) интегральная кривая системы уравнений (7) - (8) выходила из особой точки типа «седло» и приходила в точку с определенной насыщенностью (или концентрацией окислителя) при температуре зажигания Т*. Если жидкая фаза неподвижна (кокс/кероген), то рассматриваемая здесь задача относится к теории фильтрационного горения. В этом случае можно указать достаточно простые критерии реализации ИТР- или ИЬР-режимов. В обоих случаях особая точка, соответствующая максимальной температуре и исчерпанию одного из реагентов, является «седлом». Если жидкая фаза подвижна, то критерии, позволяющие различить ИТР- или ИЬР-режимы, усложняются. Для их упрощения мы наложим ограничения на исходную насыщенность жидкой фазы, которые, впрочем, не являются обременительными. Вообще говоря, уравнения (1) допускают решения, в которых насыщенность претерпевает разрывы, в то время как концентрация и температура Т непрерывны. Решения, содержащие разрывы насыщенности, следует допустить и для системы уравнений (7) - (8). Мы не будем рассматривать такие решения, а ограничимся решениями, в которых все переменные непрерывны. Тогда в уравнении (7), описывающем структуру волны горения, знаменатель всегда положителен. Так как мы считаем «г ^ в о, то на скорость фильтрации и диапазон изменения насыщенности накладывается ограничение кр'(з) < 1 при «г ^ в ^ во.
Чтобы составить качественное представление о величине скорости и структуре волны горения, рассмотрим упрощенную модель волны в ИТР-режиме. Согласно классической теории горения [11] (в частности фильтрационного) фронт волны горения может быть разделен на две зоны - зону реакции и зону прогрева. Это разделение является условным, но проявляется тем более четко, чем больше энергия активации, а точнее число Зельдовича 5(То - Тг)
/е =-2-. В зоне реакции происходит сгорание топлива и выделение тепла. Скорость
То
реакции вне этой зоны пренебрежимо мала. Вообще говоря, если волна горения распространяется нестационарно, температура в зоне реакции, которую обозначим Ть, может превышать То. Это различие Ть и То относительно невелико, но в силу сильной зависимости скорости реакции от температуры может оказать существенное влияние на скорость волны и привести к ее неустойчивости.
Выделим зону реакции по какому-либо признаку, например, по уровню насыщенности. На правой границе зоны реакции насыщенность чуть меньше 8о, на левом - чуть больше «г. Обозначим Д(£) движущийся отрезок между этими уровнями - зону реакции. Интегрируя первое уравнение системы (1) по зоне реакции, получаем
[ в(1х = (рШ - тИв)ь - (рШ -тИв)К - аг [ Кйх. (11)
М Уд Уд
Здесь {¡Ш — тБз) ь - поток нефти через левую границу зоны реакции, {¡Ш — тБз) к -поток через правую границу.
В рассматриваемом приближении принимается, что 8ь = §1, 8и = 8о и скорости левой и правой границ зоны реакции совпадают: Иь = Ии = О. Кроме того, принимается, что волна горения распространяется квазистационарно, так что левой частью (11) можно пренебречь. В итоге получаем
(цШ — тИ81) — (щШ — тОо) = аЛ К(1х. (12)
Jд
Аналогично, интегрируя по зоне реакции комбинацию первых уравнений системы (1), получаем
тИ^1 — (1 — 81) У1) — Ш (Ац — (1 — ¡1) ш) = А(тО8о — Ш<ро). (13)
К зоне реакции примыкает зона прогрева, где под действием теплопроводности температура увеличивается от значения Т до температуры Т*, при которой скорость реакции становится существенной. В зоне прогрева 8 = «1, а распределение температуры описывается уравнением
д ( \ д { дТ\
— тсТ + (1 — т) с3Т+ — ШсТ — Х— = 0. (14)
" ' у ■ дх\" ~ " дх
При этом опять принимается, что волна горения распространяется квазистационарно, т.е. зона прогрева успевает подстраиваться под изменения скорости распространения волны.
В этом приближении ширина зоны прогрева составляет = —---. На краю зоны
с( Ш — тик*)
т^т (дТ\ Т* —Т1
прогрева и в зоне реакции, где Т ^ Т, положим —— « —--.
\дх ) г ки Интегрирование уравнения энергии по зоне реакции дает
4 / Ейх = ( сТШ + — ИЕ — Х^ ) — (сТШ + — ИЕ — Х^) , (15) М Уд \ дх/ ь \ дх) и
где Е = тсТ + тдв + (1 — т)с3Т — плотность энергии.
В рамках рассматриваемого приближения считается, что Ть ~ Т* ~ Ти ~ То ^ Ть. Так как волна горения распространяется квазистационарно, левой частью (15) и потоком тепла через правую границу зоны реакции можно пренебречь. В итоге
(дТ\ т Т
~дхс) ~ ХтЬ—Т. (16)
Для вычисления интеграла по зоне реакции в правой части (12) перейдем от интегрирования по пространственной координате к интегрированию по температуре, считая, что в зоне реакции (точнее той ее части, где сгорает основной объем нефти) градиент температуры
Ть — Т1 тд
составляет —--. Имеем
I- ш
ai I Rdx & — Г fс-£/т dT - lth ai Г fe-£/TdT ai /д RdX ~ r Jt. fe (dT/dx)r =Tb -T1T k fe
kh ai fb e-£ T TbL =_^__Tb_
Tb -Ti т £ (W -k*mD)cn £(Tb -Ti)
Здесь величина ть = —— с£/Тъ характеризует время сгорания топлива в зоне реакции, Д — а1 ь
«средний» по зоне реакции предэкспонециальный множитель. Объединяя соотношения (12), (13), (16), находим
W -k*mD -\fe -Vb' (17)
т.В =С(Ть -,Г')(1 + А/"г'/" - 1+Р° -Щ. (18)
К (1 - Ро) - 1 + «о
Исключение из соотношений (17) и (18) скорости волны горения дает неявную зависимость температуры перед фронтом Ть от входных данных в стационарном режиме: темпа закачки газа Ш, концентрации в нем окислителя уг, температуры Тг и насыщенности «о перед фронтом. Время сгорания топлива в зоне реакции ~ ть. За это время фронт пройдет расстояние 1Г ~ Ить, которое можно интерпретировать как ширину зоны реакции. Ширина зоны прогрева составит = Л/(сЩ).
Скорость распространения волны горения в ИЬР-режиме оценивается аналогично и составляет
Ш -к'т° = V^ Ш-Ть? = -,Ь' (19)
или
ти = <Ы -ЪШ -ро •Щь. (20)
8о - К*ро
Здесь То — температура перед фронтом. Так же, как и выше, исключение из соотношений (19) и (20) скорости волны горения дает неявную зависимость температуры за фронтом Ть от входных данных в стационарном режиме.
Л<7 Т*
4. Режим Reaction-trailing front
Рассмотрим режим Reaction-trailing front, когда волну реакции опережает тепловая волна. В этом режиме волна горения продвигается по предварительно прогретой среде (скелет, нефть, инертные газы). Исходная насыщенность нефти в слое - sо. Концентрация кислорода в закачиваемом в слой газе равна у\, температура - Ti. Темп закачки описывается скоростью фильтрации газовой фазы Wg. В RTF-режиме температура продуктов реакции больше температуры закачиваемого в слой газа. Зададим значение температуры перед фронтом реакции То и выясним, при каких значениях параметров (темп закачки и концентрация окислителя) реализуется такой режим распространения волны горения. Для анализа RTF-режима удобно выразить константы в правых частях (6), (8) через значения переменных на правой «бесконечности». В рассматриваемой модели считается, что достигаемые при горении температуры достаточно велики, так что при таких температурах нефть и окислитель не могут сосуществовать за пределами зоны горения. Применительно к RTF-режиму это означает, что у о = 0. С учетом этого замечания уравнения (6), (8) примут вид
As - (1 - s)y - к(Ар - (1 - p)у) = А(so - kifio), (21)
Л dT
(к* - к)сТ - q(kLp - s) + —d^ = (к* - к)сТо - q(kpo - so). (22)
В RTF-режиме s1 < s0, p1 < p0, T1 < T0. Переписав (10) в виде
к - S1 - S0 (ГГ rj, S
к p1 -p0 _ C(T1 - T0)
к - к* q(p1 - P0)
< 0,
заключаем, что к лежит на отрезке, ограниченном числами к* и —-—. Согласно при* шг - шо
нятому условию ограничиться рассмотрением только непрерывных решений ( кр < 1), заключаем, что к > к*.
Соотношение (21) перепишем в виде
* = • = У-++Т- • . (23)
Р-Р У + А у + А
Из геометрических соображений ясно, что условие s ^ 0 при 0 ^ у ^ у\ будет выполнено, если
1 - sq a(yi) - SQ < ^ < a(yi) yi + ASq
(24)
1 - ро р(ух) - ро Р(уг) уг + Аро'
На рис. 2 слева показан график функции Баклея (синяя сплошная линия) на плоскости (в,(р). Точка С с координатами (з0,р0) изображает состояние жидкой фазы перед волной горения. Точки с координатами (а(у), Р(у)) лежат на отрезке (пунктирная линия), соединяющем точку С с точкой (1,1). Концентрации окислителя уо = 0 соответствует точка С, так как а(0) = Зо, Р(0) = ро. Концентрации окислителя в закачиваемом газе у\ соответствует точка В с координатами (а(уг), Р(уг)). Каждой точке на отрезке С В, параметризованном концентрацией окислителя 0 ^ у ^ у\, соответствует точка на графике функции Баклея с координатами (в,р). Согласно (23) наклоны всех хорд (красные линии), соединяющих соответствующие точки, одинаковы и равны 1/к. Если, как мы полагаем, ^ ^ Зо, то
Р(уг) - 0 уг + Асро
минимальный наклон семейства хорд составляет
P(Vi) - <Pq = 1 - <PQ
a(yi) - Sq 1 - SQ '
a(yi) - 0 yi + ASq
а максимальный
Рис. 2. (я, ^-диаграмма для анализа волны горения: ИХЕ-режим — слева, ИЬЕ-режим — справа
Таким образом, безразмерная скорость фильтрации лежит в пределах Уг + ^«о
и п!^-режим может реализоваться при выполнении условия
к± < к ^
yi + Аро
, уг + Аво тс
к* < --—, иначе говоря: уг < ---—— ро) - Аро. Отметим, что если жидкая
Уг + Аро (1 - т)с8
фаза неподвижна (р = 0), то ИТЕ-режим может реализоваться при выполнении условия тс .
Для более аккуратного расчета скорости волны горения преобразуем (21) к виду
У = А
S - Sq + к(ро - р)
1 - s - к{1 - р) :
а из (22) находим
dT_ dz
mD
(с(к - h)(Т - То) - q(s - Sq) + kq(p - Pq)^J.
(25)
(26)
Уравнения (7) и (26) образуют автономную систему ОДУ, фазовый портрет которой на плоскости (в, Т) описывается уравнением
dT (mD)2T (с(к - к*)(Т - Tq) - Q(s - sq) + kq(p - ро^(1 - кр')
ds
Хал
f (s,y) exp(-S/Т)
где у вычисляется из (25).
Выберем скорость реакции К пропорциональной насыщенностям фаз и концентрации окислителя в газовой фазе, f (в, у) = 8(1 — в)у. Уравнение (27) имеет особую точку (во, То), так как при = во согласно (25) у = 0. Линеаризуя числитель и знаменатель (27) в окрестности этой особой точки, получим
(1Т = (1 — а0 — к(1 — ро)) а(Т — Т0) + Ь(з — 80) (18 «о (1 — ) 5 — «0 '
где а = к — к*, b = 1 (Ы - 1), С = ^^ > 0, с Ало,1
то = те
£/То
Согласно сформулированным выше условиям реализации ИТР-режима, коэффициент при Т — То в числителе последнего уравнения (1 — ^ — к(1 — щ))(к — к*) < 0, так что особая точка (во, То) имеет тип «седло».
si so
Рис. 3. Фазовый портрет системы на (s, 0)-плоскости: RTF-режим — слева, RLF-режим — справа
С
Определим безразмерную температуру в = -(То — Т). Уравнение (27) перепишется в виде
лл ((к — к*)в + (s — s0) — k(p — ро Л (1 — W)
Da * = --та—1-
где *№=«*>( ^
Da =
Ха
1
> 0 — число Дамкелера, а у вычисляется из (25).
Ш 2 сто
Зафиксируем параметры среды (скелета и флюидов) т, с, с3, X, к*, реакции т, 8,д и граничные условия . Следует зафиксировать безразмерную скорость фильтрации к, но удобнее фиксировать насыщенность за фронтом волны горения «1, а связанный с ней параметр к вычислять по (9). Эти параметры должны удовлетворять условию существования бегущей волны. Находим такое значение числа Дамкелера (в размерных единицах - расхода Ха-\
Ш = \ —- ), при котором сепаратриса (красная кривая на рис. 3 слева), выпущенная
V Ба ■сто
из особой точки - «седла» (зо,0о), соединит ее с точкой («1,^1). Скорость фронта волны горения составит И = Ш/(кт).
5. Режим Reaction-leading front
Рассмотрим теперь режим, когда волна реакции опережает тепловую волну (Reaction-leading front), при этом температура за фронтом волны реакции Ti больше температуры перед фронтом Tq . Так как в этом случае содержание окислителя в зоне высокой температуры у1 > 0, считаем, что за фронтом прореагировало топливо s1 = 0.
Дальнейший анализ во многом аналогичен п. 4, но за отсчетный уровень удобно принять состояние за волной (индекс «1»). В ИЬР-режиме 0 = < Зо, 0 = р1 < ро, Т\ > То. Переписав (10) в виде
к_ 81 - 5о , ,
К Р1 - ро _ ГЛТ1 - то)
к - к* ч(р1 - Ро)
> 0,
, 1 в1 - п
заключаем, что к лежит за пределами отрезка, ограниченного числами к* и -. В
Р1 - Ро
рамках принятого ограничения к < к*. Соотношение (6) перепишем в виде
в — а , . у — у1 . .
к =-, а(у) = -—Ц- < 0. (28)
р - а у + А у ;
Из геометрических соображений ясно, что условие в ^ 0 при 0 ^ у ^ у1 будет выполнено,
, «о - а(0) у1 + Аво
если к ^--- =--—.
ро - а(0) У1 + Аро
На рис. 2 справа показан график функции Баклея (синяя сплошная линия) на плоскости (8,<р). Точка А с координатами (яо,1ро) изображает состояние жидкой фазы перед волной горения. Точка С с координатами (0, 0) изображает состояние жидкой фазы за волной горения. Точки с координатами (а(у),а(у)) лежат на луче, выходящем из точки С в третий квадрант (пунктирная линия). Концентрации окислителя перед фронтом волны у = уо соответствует точка В с координатами (а(уо),Р(уо)). Концентрации окислителя в закачиваемом газе у1 соответствует точка С с координатами (0, 0), так как а(у1) = 0. Каждой точке на отрезке С В, параметризованном концентрацией окислителя 0 ^ уо ^ у ^ У1, соответствует точка на графике функции Баклея с координатами (в,р). Согласно (28) наклоны всех хорд (красные линии), соединяющих соответствующие точки, одинаковы и равны 1/к. Если, как мы полагаем, 0 ^ уо ^ у ^ Ш, то минимальный наклон семейства хорд составляет Ро - а(0) = У1 + Аро ¿о - а(0) У1 + Аво '
В случае подвижной жидкой фазы безразмерная скорость фильтрации лежит в пре-У1 + Аво
делах --— ^ к < к* и КЬ^-режим может реализоваться при выполнении условия
У1 + Аро
, у1 + Аво тс
к* > --—, иначе говоря, у1 > ---—- ро) - Аро. Отметим, что если жидкая
У1 + Аро (1 - т)с3
фаза неподвижна (р = 0), то ИЬР-режим может реализоваться при выполнении условия
у1 > —-т—Из соотношений (3) и (5) находим
(1 - т)с3
= А(з - кр)+ У1(1 - к) У = 1 - 8 - к(1 - Р) , (29)
г!Т тВ / \
^ = ~^{с(к - К)(Т - П) + д(кр - з)). (30)
На фазовой плоскости (в, Т) структура фронта волны горения описывается уравнением
йТ _ (тИ)2т (с(к - к*)(Т - Т1) + ч(к(Р - 8)) (1 - V)
йв Ха1 / (з,у)ехр(-£/Т)
(31)
где у вычисляется из (29).
Особой точкой уравнения является точка (,в1,Т1). Линеаризуя числитель и знаменатель (31) в окрестности этой точки, получаем
дТ = са(Т - Т{) + Ъз йв в ,
где а = к — к*, Ь = —, С = —-- > 0, п = те£/Тг. В силу предыдущих рассужде-
с \yiai
ний особая точка (в 1,Т\) уравнения (32) имеет тип «седло» и интересующее нас решение изображается на фазовой плоскости одной из сепаратрис.
Определим безразмерную температуру в = -(Т\ — Т). Уравнение (32) перепишется в виде ^
м ((к — к*)в — к^ + 5)(1 — у) =--т
/ 8в/Т1 \ Ха1
где Ф(6М = ехр -;-- , Ба = —~— > 0 - число Дамкелера.
\сТ1/д — 0) Ш 2ст1
Аналогично п. 4 зафиксируем параметры среды т,с, с3,\, к*, реакции т, 8,д и граничные условия во,у1. Следует зафиксировать безразмерную скорость фильтрации к, но удобнее фиксировать концентрацию окислителя перед фронтом волны горения , а связанный с ней параметр к вычислять по (9). Эти параметры должны удовлетворять условию существования бегущей волны. Находим такое значение числа Дамкелера, при котором сепаратриса (красная кривая на рис. 3 справа), выпущенная из особой точки — «седла» (в 1, $1), соединит ее с точкой (во,&о). Скорость фронта волны горения составит V = Ш/(кт).
6. Численное моделирование распространения волны внутрипластового горения
Для численного моделирования распространения волн внутрипластового горения был разработан вычислительный алгоритм, основанный на расщеплении по физическим процессам. На первом этапе рассчитывалось распределение давления, насыщенности и концентрации окислителя. Для этого численно интегрировались уравнения баланса (1) с помощью неявной противопоточной консервативной схемы. На втором этапе интегрировалось уравнение энергии (2) по явной (по температуре) схеме.
Рассматривался пористый слой длины Ь = 10 м. Шаг сетки по пространству выбирался из условия, что зона реакции, ширина которой оценивалась по формулам п. 3, должна содержать примерно 10 расчетных точек, и составил Ъ = 510-3 м. Шаг по времени составил 0.3 с. Была выбрана квадратичная зависимость относительных фазовых проницаемостей от соответствующих насыщенностей, при которой функция Баклея имеет вид
= кь (з)/¡¡ь = в2 кь(в)/ль + кс(«)/лс 82 +е(1 — в)2,
где ¡ь, ¡с — коэффициенты динамической вязкости жидкой и газовой фаз, е = ¡ь/лс.
Далее опишем моделирование распространения волны горения в ИТР-режиме.
Начальные условия на интервале 0 ^ х ^ Ь были выбраны в виде: у(х) = 0, в(х) = «о = 0.25, Т(х) = 375 К. На входной (левой) границе слоя фиксировались значения скорости фильтрации Ш = Шс = 2.63 ■ 10-5 м/с, концентрация окислителя (в газовой фазе) у = у1 = 0.231, насыщенность 8 = 0, температура Т = Т! = 367 К. На выходной (правой) границе слоя задавалось давление. Так как все теплофизические характеристики флюидов и скелета не зависят от давления, его конкретная величина не имеет значения.
Параметры задачи: стехиометрические коэффициенты а1 = 1, а2 = 10, а3 = 11, пористость т = 0.3, теплоемкость жидкой и газовой фаз с = 2.49 ■ 107 Дж/(м3-К), теплоемкость скелета с3 = 4.35 ■ 107 Дж/(м3-К), калорийность нефти д = 3.74 ■ 1010 Дж/м3, коэффициент теплопроводности А = 3 Вт/(м-К), температура активации 8 = 10 000 К, характерное время реакции т = 2 ■ 10-7 с, /( 8, у) = «(1 — 8 )у, при Т < 369 К скорость реакции К полагалась равной нулю. Отношение коэффициентов вязкости фаз е = 20. При этих параметрах к* = 5.08, А = 10, (у1 + Аз0)/(у1 + А^р0) = 9.54 > к*, так что может реализоваться ИТР-режим.
Согласно решению вида «бегущая волна» при указанных параметрах среды и «граничных условиях», температура перед фронтом составила То = 485 К, значение насыщенности
жидкой фазы за фронтом - 51 = 0.05, число Зельдовича - 2е = 5, ширина зоны реакции -1Г ~ 15 см., скорость фронта волны горения И составила 1.231 • 10-5 м/с. Величина предэкс-поненциального множителя з(1 — в)у не превышает 0.25 • (1 — 0.25) • 0.231 = 0.0433, так что скорость волны, вычисленная по приближенной формуле (18), не превышает 1.38940-5 м/с. Безразмерная скорость фильтрации к = Ш/(тБ) = 7.122.
Результаты численного моделирования представлены на рис. 4, 6 (в момент £ = 100 ч). На рис. 4 показано распределение температуры, на рис. 6 — насыщенности (черным) и концентрации окислителя (красным). Впереди распространяется тепловая волна, за ней — волна горения. Скорость волны горения в вычислительном эксперименте составила 1.234 • 10-5 м/с.
Рис. 4. Распределение температуры Т(ж) в ИТЕ-режиме
Рис. 5. Распределение температуры Т(ж) в ИЪЕ-режиме
Рис. 6. Распределение насыщенности жидкой фазы «(ж) и концентрации кислорода у(х) в ИТЕ-режиме
Рис. 7. Распределение насыщенности жидкой фазы в (ж) и концентрации кислорода у(х) в ИЪЕ-режиме
При моделировании распространения волны горения в КЬР-режиме на входной (левой) границе слоя фиксировались значения скорости фильтрации Ш = Шс = 3.1740-5 м/с, концентрация окислителя (в газовой фазе) у = у1 = 0.231, насыщенность 5 = 0, температура Т = 380 К. На выходной (правой) границе слоя задавалось давление.
Начальные условия, параметры среды и кинетики не отличались от предыдущего случая за исключением того, что теплоемкость скелета составляла с8 = 1.3 • 108 Дж/(м3•К), а скорость реакции Я полагалась равной нулю при Т < 379 К. При этом к* = 13.12, А = 10, (у1 + )/(у1 + ^^о) = 9.54 < к*, так что может реализоваться ИЬР-режим.
Согласно решению вида «бегущая волна» при указанных параметрах среды и «граничных условиях», температура перед фронтом равнялась начальной температуре в слое То = 375 К, температура за фронтом составила = 485 К, значение концентрации окислителя (в газовой фазе) перед фронтом — у 1 = 0.05, число Зельдовича — 2е = 4.68, ширина зоны реакции — 1Г « 13 см, скорость фронта волны горения И составила 1.008 • 10-5 м/с. Величина предэкспоненциального множителя з(1 — в)у не превышает 0.0433, так что скорость волны, вычисленная по приближенной формуле (20), не превышает 1.17 • 10-5 м/с. Безразмерная скорость фильтрации к = Ш/(тБ) = 10.48.
Результаты численного моделирования представлены на рис. 5, 7 (в момент I = 220 ч). Впереди распространяется волна горения, за ней — тепловая волна. Скорость волны горения в вычислительном эксперименте составила 1.015 • 10-5 м/с.
7. Заключение
В данной работе рассмотрены два режима распространения волны горения при двухфазном течении в пористом слое. Установлены критерии, необходимые для реализации этих режимов. Проведено численное моделирование распространения волн горения в обоих режимах, продемонстрировавшее справедливость предложенных критериев.
Работа выполнена в лаборатории флюидодинамики и сейсмоакустики МФТИ при финансовой поддержке Российского научного фонда (проект 15-11-00015). Авторы благодарят
А. В. Колдобу, А. В. Конюхова и Ю.И. Скалько за полезные обсуждения.
Литература
1. Алекперов В.Ю., Грайфер В.И., Николаев Н.М. [и др.] Новый отечественный способ разработки Баженовской свиты (часть 1) // Нефтяное хозяйство. 2013. № 12. C. 100105.
2. Алекперов В.Ю., Грайфер В.И., Николаев Н.М. [и др.] Новый отечественный способ разработки Баженовской свиты (часть 2) // Нефтяное хозяйство. 2014. № 1. C. 50-53.
3. Добрего К.В., Жданок С.А. Физика фильтрационного горения газов. Минск: Институт тепло- и массообмена им. А.В. Лыкова НАНБ, 2002.
4. Kristensen M.R., Gerritsen M.G., Thomsen P.G. [et. al.] An Equation-of-State Compositional In-Situ Combustion Model: A Study of Phase Behavior Sensitivity // Journal of Porous Media. 2009. V. 76, I. 2. P. 219-246.
5. Souza A., Marchesin D., Akkutlu I. Wave sequences for solid fuel adiabatic in-situ combustion in porous media // Comp. Appl. Math. 2006. V. 25, I. 1. P. 27-54.
6. Mailybaev A.A., Bruining J., Marchesin D. Analysis of in situ combustion of oil with pyrolysis and vaporization // Combustion and Flame. 2011. V. 158, I. 6. P. 1097-1108.
7. Youtsos M.S.K., Mastorakos E. Numerical simulation of thermal and reaction waves for in situ combustion in hydrocarbon reservoirs // Fuel. 2013. V. 108. P. 780-792.
8. Gargar N.K., Mailybaev A.A., Marchesin D., Bruining J. Compositional effects in light/medium oil recovery by air injection: vaporization vs. combustion // Journal of Porous Media. 2014. V. 17, I. 11. P. 937-952.
9. Завьялов И.Н., Конюхов А.В. Экспериментальное и численное исследование неустойчивости многофазной фильтрации при наличии реакции с образованием газовой фазы // Материалы международной конференции «Нелинейные задачи теории гидродинамической устойчивости и турбулентность». 2014. С. 89-91.
10. Бурже Ж., Сурио П., Комбарну М. Термические методы повышения нефтеотдачи пластов. М.: Недра, 1988.
11. Зельдович Я.Б., Баренблатт Г.И., Либрович В.Б., Махвиладзе Г.М. Математическая теория горения и взрыва. М.: Наука, 1980.
12. Желтов Ю.П. Механика нефтеносного пласта. М.: Недра, 1975.
References
1. Alekperov V.Yu., Grayfer V.I., Nikolaev N.M. [et. al.] New Russian oil-recovery method for exploiting the Bazhenov Formation's deposits (part 1). Oil industry. 2013. N 12. P. 100-105. (in Russian).
2. Alekperov V.Yu., Grayfer V.l., Nikolaev N.M. [et. al.] New Russian oil-recovery method for exploiting the Bazhenov Formation's deposits (part 2). Oil industry. 2014. N 1. P. 50-53. (in Russian).
3. Dobrego K.V., Zhdanok S.A. Physics of filtration combustion of gases. Minsk: A.V. Luikov Heat and Mass Transfer Institute of NAS of Belarus, 2002. (in Russian).
4. Kristensen M.R., Gerritsen M.G., Thomsen P.G. [et. al.] An Equation-of-State Compositional In-Situ Combustion Model: A Study of Phase Behavior Sensitivity. Journal of Porous Media. 2009. V. 76, I. 2. P. 219-246.
5. Souza A., Marchesin D., Akkutlu I. Wave sequences for solid fuel adiabatic in-situ combustion in porous media. Comp. Appl. Math. 2006. V. 25, I. 1. P. 27-54.
6. Mailybaev A.A., Bruining J., Marchesin D. Analysis of in situ combustion of oil with pyrolysis and vaporization. Combustion and Flame. 2011. V. 158, I. 6. P. 1097-1108.
7. Youtsos M.S.K., Mastorakos E. Numerical simulation of thermal and reaction waves for in situ combustion in hydrocarbon reservoirs. Fuel. 2013. V. 108. P. 780-792.
8. Gargar N.K., Mailybaev A.A., Marchesin D., Bruining J. Compositional effects in light/medium oil recovery by air injection: vaporization vs. combustion. Journal of Porous Media. 2014. V. 17, I. 11. P. 937-952.
9. Zavyalov I.N., Konyukhov A.V. Experimental and numerical study of unstable multiphase filtration which occurs with gas phase creation. Proceedings of the international conference «Nonlinear problems in the theory of hydrodinamic stability and turbulence». 2014. P. 8991. (in Russian).
10. Burzhe Zh., Surio P., Kombarnu M. Thermal methods of enhanced oil. М.: Nedra, 1988. (in Russian).
11. Zeldovich Ya.B., Barenblatt G.I., Librovich V.B., Makhviladze G.M. The mathematical theory of combustion and explosion. Moscow: Nauka, 1980. (in Russian).
12. Zheltov Yu.P. Mechanics of an oil reservoir. Moscow: Nedra, 1975. (in Russian).
Поступила в редакцию 04.08.2016
