CC BY
7
Вычислительные технологии
Том 23, № 6, 2018
Расщепление по физическим процессам в некоторых задачах флюидодинамики в гидратосодержащих пористых средах
П. И. Рлгимли1, Ю. С. ШАРОВА2'*, О. Р. Рлгимли1, В. О. ПодрыгА2,3, И. В. Гасилова2, С. Б. Попов2, Ю.А. ПОВЕЩЕНКО2'4
1 Московский физико-технический институт, Долгопрудный, Московская область, Россия 2Институт прикладной математики им. М.В. Келдыша РАН, Москва, Россия 3Московский автомобильно-дорожный государственный технический университет, Москва, Россия
4Национальный исследовательский ядерный университет "МИФИ", Москва, Россия *Контактный e-mail: [email protected]
Для развития технологий освоения газогидратных залежей необходимы исследование и моделирование проблем фазовых превращений газогидратных включений в пористых средах. В данной работе с помощью предложенной двублочной математической модели диссоциации газовых гидратов в пористой среде, основанной на расщеплении по физическим процессам, проведены расчеты водонасыщенности, растепленности, а также термодинамических параметров (давления и температуры). Численная модель позволяет дискретизировать задачу в одномерном случае и реализовать абсолютно устойчивую разностную схему. Полученные результаты демонстрируют применимость предложенной модели для решения типичных задач газогидратной флюидодинамики, в том числе для исследования сложной динамики водо- и гидратонасыщенности пласта с учетом адиабатического расширения газа в коллекторном пространстве.
Ключевые слова: газовые гидраты, фильтрация, математическое моделирование.
Библиографическая ссылка: Рагимли П.И., Шарова Ю.С., Рагимли О.Р., По-дрыга В.О., Гасилова И.В., Попов С.Б., Повещенко Ю.А. Расщепление по физическим процессам в некоторых задачах флюидодинамики в гидратосодержа-щих пористых средах // Вычислительные технологии. 2018. Т. 23, № 6. С. 64-79. БОТ: 10.25743/1СТ.2018.23.6.007.
Введение
Наличие водяных паров в углеводородных газах с повышенным давлением при низких температурах приводит к их конденсации и образованию газосодержащих ледяных пробок (газовых гидратов), что осложняет процессы как транспортировки, так и переработки углеводородов. Исследования показали, что образование газовых гидратов
© ИВТ СО РАН, 2018
также возможно в природных условиях, причем объемы природных залежей газовых гидратов крайне велики [1-4]. Возможность высвобождения углеводородов из газовых гидратов делает данные соединения крайне перспективным объектом для исследования в качестве альтернативного источника углеводородов. Для этого нужны детальные исследования, направленные на развитие технологий добычи и решения возможных экологических проблем, в том числе связанных с выбросами газа в атмосферу.
Для исследования фильтрации с учетом диссоциации газовых гидратов могут использоваться уравнения механики сплошной среды, выражающие законы сохранения массы, импульса и энергии [5]. Математическое моделирование подземных газовых гидратов имеет давнюю историю и осуществляется в разных направлениях. В ряде работ в одномерном случае путем введения упрощающих предположений задача приводится к автомодельному виду. В результате выбора соответствующих переменных система уравнений в частных производных преобразуется в систему обыкновенных дифференциальных уравнений и исследуется аналитически и численно.
Подобный подход использовался [6-15] для решения ряда задач. В работе [6] рассматривается схема подвода тепла закачкой горячей воды под подошву гидратного пласта. В [7] на основе классической задачи Стефана были предложены математические модели разложения гидратов в пористой среде при депрессионном и тепловом воздействии и получены точные балансовые уравнения на фронте разложения. Фильтрация газа рассматривалась в изотермической и баротропной постановке, изменение равновесных условий гидратообразования в процессе истощения не учитывалось. В работе [8] решается задача разложения гидрата аналогично классической задаче Стефана с резким фронтом разложения.
Математическая модель с объемным характером разложения гидратов получила дальнейшее развитие в [10]. Результаты математического моделирования течения в пористых средах с фазовыми превращениями в месторождениях природного газа, содержащего газовые гидраты изложены в [11]. В работе [12] приведен алгоритм решения задачи неизотермической фильтрации пластовых флюидов в пористой среде с учетом фазовых превращений газовых гидратов, позволяющий установить степень влияния кондуктивного и конвективного переноса тепла, а также эффекта Джоуля — Томсона на характер распределения искомых параметров процесса (давления, температуры, гид-ратонасыщенности) при различных начальных и граничных условиях.
В работе [13] для оценки интенсивности оттаивания и сопутствующего газовыделения при добыче нефти и газа на северных месторождениях предложена модель теплового взаимодействия добывающей скважины и толщи многолетнемерзлых пород, содержащей реликтовые метастабильные гидраты. Получены автомодельное решение задачи и аналитическая зависимость для определения границы фазового перехода, выведена формула распространения радиуса теплового воздействия скважины. В рамках автомодельного решения проанализировано влияние гидратонасыщенности породы и теплопроводности цементного камня на радиус оттаивания грунта скважины.
В работе [14] представлены результаты теоретического исследования процесса ин-жекции газообразного диоксида углерода в пористую среду, насыщенную в исходном состоянии метаном и льдом, предложена математическая модель тепломассопереноса в пористой среде, сопровождающегося образованием гидрата диоксида углерода. Получены решения, описывающие распределение температуры и давления в пласте. Проанализированы условия, при которых реализуются различные режимы процесса образования гидрата. Построены диаграммы существования данных режимов.
В работе [15] проведено исследование инжекции углекислоты в жидкой фазе в истощенное месторождение природного газа. Предложена математическая модель процесса, учитывающая образование гидрата CO2 и вытеснение метана. Найдено асимптотическое решение задачи в одномерном приближении. Показано, что для термодинамически непротиворечивой постановки задачи от предположения, что диссоциация гидрата происходит в узкой зоне подвижной границы, следует перейти к рассмотрению области объемных фазовых переходов [9]. В этой работе задача разложения гидратов в пористой среде изучалась на основе совместного решения уравнения фильтрации газа и теплопе-реноса. В результате автомодельного решения системы получены распределения давления и температуры в пласте. Однако решения выполнялись без учета движения воды и влияния газа на изменение температуры. Аналогичного подхода в исследованиях декомпозиции гидрата в пористой среде придерживались и зарубежные авторы [16-18].
Основу кинетических моделей диссоциации гидрата составляет уравнение Кима и Бишной [19], связывающее количество выделившегося газа из гидратов с изменением термодинамических параметров — давления и температуры. Задача о разложении гидрата в этом случае сводится к системе уравнений с дополнительными источниками массы в правой части [20]. Так как численные методы решения такой системы хорошо известны [21], эта модель получила широкое распространение [22-24].
Для сравнения различных флюидодинамических моделей с газогидратными включениями на базе The National Energy Technology Laboratory и The U.S. Geological Survey (США) проводятся международные исследования [25]. Постоянно возникают новые методы, например, в Германии — SUGAR на основе PetroMod [26], в Норвегии — Retraso-CodeBright (RCB) [27]. Но в опубликованных работах не представлена методика совместного решения систем уравнений, описывающих процессы в различных областях, каждая из которых характеризуется собственным набором сосуществующих фаз, а согласование вычислительных схем для них не является автоматическим процессом. Поэтому разработка отечественного математического и программного обеспечения для решения подобных задач является актуальной задачей.
Исследование автомодельных решений позволяет подробно и глубоко изучить некоторые модельные задачи, которые во многих случаях имеют прямое практическое значение (например, для анализа работы скважин), но для расчета более сложных задач необходимо использование численных методов. Это важно также при исследовании одномерных задач в более сложной постановке. Такие расчеты проводились в [28, 29] и др., но методы, предлагаемые в данной работе, позволяют расщепить перенос процессов насыщения и диссипативный блок и использовать неявную схему только для давления. Кроме того, они дают возможность проводить единый расчет во всей плоскости Р, Т.
В настоящей работе в качестве математической модели используется система уравнений [30-32], в которой наиболее полно учтены основные физические особенности процесса квазиравновесной многофазной фильтрации при наличии газогидратных включений. Данная система уравнений сводится к блоку переноса процессов насыщения, отвечающему за конвективный перенос насыщенностей Sw, Su и обладающему в основном гиперболическими свойствами и дополняющему этот блок диссипативному уравнению пьезопроводности для давления. Такой подход позволяет применять явно-неявные разностные схемы при решении задачи (неявные только для давления). Разработан численный метод решения пространственно-одномерной задачи, базирующейся на разработанной математической модели, проведен ряд численных экспериментов, в том числе с учетом адиабатического расширения газа в коллекторе.
1. Описание математической модели
Математическая модель фильтрации состоит из блока переноса процесса насыщения и диссипативного уравнения пьезопроводности, записанного для давления. Блок переноса процессов насыщения имеет вид
д_ дЪ
— {га(ЯРп + (1 - я)РиРт)} + [РтV™] + = 0,
д_ т
Диссипативное уравнение пьезопроводности записывается как
— {т(Б„ (1 - Бад )рд + (1 - Я )Ри (1 - Рт))} + ^У [Рд Vg ] + ^ =0.
дР ф
от тр^
Физика уравнения (3) заключается в материальных коэффициентах:
(Ри )р . тр
Бр = т8А Б,
к^ + (1 - ^) + (1 - + ЩЛ +
I [. Рад Рд \ Р» т J
Ф
+--{т[8и Бад р ад (^ад )р + Б„ (1 - Бад )Рд (£д )р + (1 - Б„ )р„ (£ у )р] + [(1 - Ш) Р 3£ 3]})}
Рад^ад + (1 Рад ^д ,
Ф /Зад + (1 - /Зад) 1
три Рад Рд ри
В1С = — [Рад V,« ] + - а1у [рд Vg ] + ^ + ^,
Рад Рд Рад Рд
БЮе = РадУад ▽ £ад + Рд Vg ▽ £д + ^У [Р (V ад + Vg)] + diУ W + (^ £ &ад Чад
W = -{т[Б„(БадХад + (1 - вад)Хд + (1 - Б„)X,] + (1 - т^} ▽ Т. Здесь Ир — барический коэффициент; 6£ — скачок внутренней энергии среды при фазо-
Ф
вом переходе, отнесенный к единице массы;--удельный скачок объема. Индексы
три
д, 'ш, и, в относятся к газу, воде, гидрату, скелету пористой среды; I — индекс, указывающий фазу; Р — давление; Т — температура; Бад — водонасыщенность; Бд = 1 - Бад — газонасыщенность; V — гидратонасыщенность; Би = 1 - и — растепленность; р\ = (Р, Т) — плотность фаз; VI — скорость фильтрации соответствующей фазы; /Зад — массовая доля воды в гидрате; т(г, Р) — пористость; г — радиус-вектор; £ — время; ф (Ь, г, Бад,Б„, Р) — плотности источников фаз; £г(Р,Т) — внутренние энергии фаз; Хг(Р,Т) — коэффициенты теплопроводности, входящие в тепловой поток W.
Скорости фильтрации даются выражениями (закон Дарси с учетом гравитации в среде с общим давлением):
^^ад = - (▽Р - драд к), Vg = - ^^ (▽Р - дрд к),
Рад Рд
где к — вертикальный координатный орт; д — ускорение свободного падения; к(г, Би, Р) — абсолютная проницаемость; кгад(Бад), кгд(Бад) — фазовые проницаемости; рад, рд — вязкости воды и газа.
Внутренняя энергия гидрата выражается через энергии создающих его газа и воды:
Р'ш+ (1 Р'ш + К,
Р
где К — скрытая теплота фазового перехода единицы массы гидрата, ц = £1 +---
Р1
энтальпия.
Состояние гидрата описывается соотношением фазового равновесия
Т = А 1п Р + В, (4)
где А, В — эмпирические константы. В силу этого соотношения зависимость от Т во всех выражениях можно свести к зависимости от Р. Уравнение состояния для газа имеет вид
_ РМ Рд = г (Р,Т )ПТ,
где М — молярная масса газа, К — универсальная газовая постоянная, г(Р, Т) — коэффициент сверхсжимаемости газа.
Таким образом, для трех уравнений (1)-(3) модели с учетом (4) неизвестными являются три независимых переменных: влагонасыщенность , растепленность и давление Р. Все уравнения (1)-(3) нелинейны.
Получившаяся система состоит из функционального блока (1), (2), отвечающего за характеристический перенос насыщенностей (в математическом плане - это гиперболичность в независимых переменных , на фоне фиксированного давления Р), и функционального блока (3), описывающего диссипативные и конвективные процессы, выраженные нестационарностью по времени первого порядка (д/дЬ) и пространственными дифференциальными операциями второго порядка (в терминах вектора ▽). В последнем случае независимой переменной является давление Р при фиксированных сатурациях и . Отметим, что перенос насыщенностей — растепленности и влагонасыщенности — вдоль характеристик связан соответственно со сносом вниз и сносом вверх по потоку.
Дивергентные слагаемые, входящие в выражение для ИЮ, имеют тот же смысл и структуру, что и в обычном уравнении пьезопроводности для двухфазной среды [5]. Однако величина ИЮ, отвечающая за конвективный перенос масс, воды и газа, а также источники воды и газа, в уравнении трехфазной пьезопроводности (3) домножается на удельный скачок внутренней энергии среды при фазовом превращении 6£. Аналогично выражение ОЮ£ домножается на удельный скачок объема при фазовом превращении ф/трг,. Величина ОЮе связана с конвективным переносом внутренних энергий воды и газа, работой сил давления, а также с диссипацией тепловой энергии и тепловыми источниками. Значение барического коэффициента Ир может интерпретироваться аналогично коэффициенту упругоемкости пласта, определяющему скорость распространения возмущений давления в пласте. Физическое значение слагаемых, входящих в Ир, уточняется в конкретных задачах.
2. Постановка задачи
Рассматривается одномерная постановка задачи. Предполагается, что доля теплопроводности в общем балансе переноса тепла пренебрежимо мала по сравнению с конвекцией, т. е. в уравнении энергии кондуктивная составляющая полагается равной нулю
(&у W = 0). Ускорение свободного падения также не учитывается (д = 0). В итоге система (1)-(3) сводится к системе следующего вида:
д д — {m(SVSWpw + (1 - Sv)pvfiw)} + —
— {m(Sv(1 - Sww)pg + (1 - Sv)pv(1 - pw))} + тт-
pwVw
+ qw = 0,
д
дх
pgVg
+ qg = 0
Ф
D„ ддр +
(
L д
pw дх
1 д
w Vw g д х g Vg
, r . лг д^-ш . Л, дед д +--pwVw+ pgVg— + —
m pv\ дх дх дх
k krw I дР \ тк krg I дР
Р (К, + Vg)
+ Ч~ + ^ +
w pg
+ (qe - ZwQw - £:
Vw =
^ V = ~rg ,
дх Г g ug \ дх
\ ОХ ) ^д
Границы расчетной области предполагаются непроницаемыми твердыми "стенками" т. е. поток через них нулевой:
Vw
х=0
0, Vg
g I х=0
0,
Vw \х=1
0,
Vg \x=i
0, t> 0.
Важным обстоятельством является то, что исходная задача, сформулированная в виде законов сохранения (массы Н20, СН4 и полной энергии среды) с общей матрицей системы относительно функций , Sw, давления Р и температуры Т, обладает смешанными гиперболическими и параболическими свойствами. Непосредственное использование такой системы для целей определения динамики переменных Su, Sw ,Р ,Т и построения неявной разностной схемы, требуемой для расчетов параболических уравнений с крупными шагами по времени, затруднительно.
В разработанном численном методе происходит расщепление на блок переноса на-сыщенностей флюидов на фоне заданного поля скоростей (обладающий в основном гиперболическими свойствами) и пьезопроводный блок системы с гидратными включениями, определяющий диссипативную эволюцию термодинамических параметров равновесной флюидогидратной модели. Такое расщепление по физическим процессам позволяет создавать прикладные алгоритмы с матрицей системы разностных уравнений только для давления Р, что дает возможность производить расчеты с меньшим числом неизвестных.
0
3. Описание вычислительного эксперимента
При решении поставленной начально-краевой задачи используется метод конечных разностей. Для этого строится сетка с пространственным равномерным шагом h и шагом по времени тп, где п — номер шага по времени, а исходные уравнения, граничные и начальные условия заменяются их сеточными аналогами. При построении схем используется аппроксимация UPWIND для влагонасыщенности и аппроксимация DOWNWIND для растепленности. Данный прием следует из анализа гиперболичности системы уравнений относительно Su, Sw на фоне фиксированного поля скоростей, определяемого законом Дарси. Для простоты верхние индексы п +1 опущены.
Итоговые разностные уравнения имеют вид
Ор(Рг) + 6£(Рг)ОЮг +(
т \rnpv)
тСег = 0, г = 1,..., N - 1, п> 0, (5)
ОЮг = - —
Рад
1
+-
Рад
к(Б^,г+1)кГад ($ад,г)ег + к(Б^,г)кГад ($ад,г+1 )(1 - ег) к(Бп,г)кГад ($адл- 1)/г + к(Б^г-1 )к-ад таа - /г)
Рг+1 - Рг
Рд,г+0.5
Рд Рд,г
Рд,г-0.5
Рд Рд,г
к($и,г+1)кгд ($ад,г)ег + к(Бп,г)кгд № )(1 - ег)
Рг - Рг-1
К2
Рг - Рг-1
К2
Рг+1 - Рг
+
)!г + к(БПг (зад,,г)(1 -1г)
К2
К2
+ ^ + ^, (6) Рад Рд,г
хА
хА
Рг+1 - Рг
К .
Рг - Рг-1
К
(
Рад Сад
РадРг+0.5
к(Бп,г)кГад ($ад,г+1)
I г + -
(
РадС-ад
РадРг-0.5
) -ад ( Б ад г )
Рд,г+0.5 Сд РдРг+0.5
Рд,г-0.5 Сд
X
к^к-д ($ад,г+1)
) -д( ^ад ,г)
X
(1 - Iг) - Рг+0.5
п
и,г+1
РдРг-0.5
)к-ад( Бадл)ег + к(Б^Кад( $адл+1)(1 - ег) )Р—1 +
+(к(БП г+-\)кгд (Бад,г) ег + к(Б1пЛ) кгд
г+1 )(1 - Сг))р
1
Рг+1 - Рг К2
+ Рг-
г-0.5
п
у,г>
)к-ад(Бадл_ 1) ¡г + к(БПг к-ад( ^¿Х1 !г) )Р— +
г=
+ (к(БПгг)к -д( $ад,г-1 ) /г + к(Б1г -1) "'-д
+ (0_е,г 0_ад,г&ад,г Яд.г^д.г^ 1, -(Рг+1 -Рг) > 0, (1, -(Рг+1 - Рг) > 0, ,
Рг- Рг
г1
К2
(7)
0, -( Рг+1 -Рг) < 0,
г
I
0, -( Рг+1 -Рг) < 0,
!
1, Рг+1 > Рг > Pг-l, 0, Рг+1 <Рг < Рг-1,
Ер,г — ш 8е,гБи,г(1 Б,ш,г)-
(Рд,г)Р , ( Ф
+ ( ) (1 -т) Рз (£ г )р+
\трV) г
.(А.
+ 1 ) {т[Би,гБад,гРад(£ад,г)р + Б„,г(1 Бад,г)Рд,г(£д,г)р + (1 $и,г)Ри(^и,г)р]}
V тР» / г
л =яг +(1-я )р -г = А ^И^Рад)-!
Ое,г Рад&ад,г + (1 Рад)£д,г &и,г, + :
^•Ри Рад Рд,г Ри
2
2
трг
С С _ СП СП 1 _ с _ /1 _ СП N
,--+ тр„(^--+
Рь
Кцъ
k(Sv,í+1)krw{^Ь,^ тах (
Т
/ туп _ туп \
(- ^ ^ °) +
+к(Кь( SW,г+l)mi^ -Рь
туп _ туп
Рг+1 о
К
К/,
рп _ рп
к(^г)кГь (3^_1)тах( - * К 1-1, 0 1 +
/ рп _ рп
+к( ^)кгь( 5Ь,<)тт - * К г_1 , о г=1,...М- 1, п> 0,
+ = 0,
+
(1 Sw,í)рд,% ^Л(1 , /-. о N 1 ^^ (1 ,
т--+ тр (1 - ^3w)--+
',¿+0.5
к(^'™г+1)кгд(^'г«,г) тах (
К/Лд
+к(S™г)кrд( -
рга _ рга
Рг+1 р
К
туп _ туп
Рг+1 0
К
0+
<Л—0.5
К/д
/ рп _ рп \
ВД,г)кГд^,г_1)тах - * К , 0^ +
/ рп _ рп
+к(З^к-д( SW,г) mi^ - * К г_1 , 0 г=1,...М- 1, п> 0,
+ д,
:ю)
р0
уп
р(Х, 0), г=1,М - 1, SWА = Sw(Х, 0), = SV(Х, 0), г = 0,М,
= 0, УдПо = 0, = 0, Удпм = 0, п> 0 (ро = рърм = рм_1).
Здесь , сд — теплоемкость воды и газа при постоянном объеме. Интерполяционные множители &±о.5 = 0.5(д^ + ^±1), где д = р, Рд, ен, ^, ^ € [0, 1].
Очевидно, что разностная схема (5)-(10) обладает первым порядком аппроксимации по пространству и времени на гладких решениях. Метод решения заключается в следующем. Сначала решается диссипативный блок (5)-(8) с фиксированными значениями растепленности и водонасыщенности. Соответствующее разностное уравнение представляет собой систему нелинейных алгебраических уравнений, которая после применения метода хорд в индексной форме сводится к трехточечному уравнению
Аг5рг_1 - Сг5рг + Вг5рг+1 = -Ъ, г = 1, М,
где 5 — номер итерации, 8рг = р/+1 - р/ — приращение давления на 5 + 1-й итерации.
Уравнение (11) на каждом временном слое решается с помощью метода простых итераций в сочетании с алгоритмом прогонки. Итерационный процесс продолжается до тех пор, пока не будет достигнута заданная точность по давлению:
< £11рП + £2,
0
где £1, £2 — малые величины.
Далее, располагая значениями давления Р, можно совершить переход ко второму этапу, где выполняются расчеты Б* и по явной схеме. Для этого из дискретных соотношений (9) и (10) после сложения находятся значения растепленности, которые затем подставляются в уравнение (9), откуда определяются значения водонасыщенно-сти. Процедура повторяется до заданного времени расчета. Машинное представление чисел удается проводить в размерном виде.
4. Задание начальных данных
Рассмотрим процесс, в начальный момент времени которого давление распределено по закону
Р(х, 0) = 2 ■ 107 Па, а водонасыщенность и растепленность также однородны по пространству:
(х, 0) = 5***,, (х, 0) = Sv,
где 0 < Б"* < 1 и 0 < Б1* < 1 — постоянные величины, х Е [0, /], 1 — длина расчетной области.
Рассмотрены два случая работы источника в среде коллектора:
^ = а(Р -Р*), дд = /3(Р -Р*), д£ = е^ + £ддд, (12)
а также источники энергии с учетом адиабатического расширения газа
Яе = + ВдЯд + ~Р, (13)
д
где Р* — постоянная величина; а и 3 — постоянные параметры, являющиеся некоторыми характеристиками источника (стока) в призабойной зоне, обусловленного перепадом давления внутри скважины и в пласте. Выбраны следующие значения параметров, характерные для Мессояхского газогидратного месторождения:
кг кг кг
—, pv = 910 —, ps = 2800 —
м3 м3 м3
Дж ОГАА Дж
pw = 103^, р„ = 910 ^г, ps = 2800 ^г, ßw = 0.9, А = 7.28 К, В =169.7 К,
3 3 3
pw = 10-3 Па ■ с, pq = 0.014 ■ 10-3 Па ■ с, cw = 4165 , с„ = 2500 ^^,
кг ■ K кг ■ K
Дж кг Дж Дж
s = 873 , М = 0.016 -, h = 514810—, R = 8.31 ^ , т = 0.35,
кг ■ К моль кг моль ■ К
SW = 0.6, S* = 0.5, k(Sv) = ko(Sv)3, ко = 10 мД = 10-14 м2, Р* = 2 МПа, krw(Sw) = 1.477SW - 1.587SW + 1.11SW - 0.0473, krg(Sw) = 1.044 - 1.7Sw + 0.6SW.
Минимальное значение водонасыщенности Sw min = °.55, krw (Sw) = 0 krg ( Sw ) = krg (Sw min) при Sw ^ Sw min. Максимальное значение в°донасЫщенности Sw max — 0.9,
krw (Sw) krw (Sw max)? krg (Sw) 0 при Sw ^ Sw max.
Длина модельной трубы полагалась равной I = 1 м, шаг по пространственной координате h = 0.01 м. Расчеты проводились для моментов времени t = 1,10, 50 с. Границы области действия источников х Е [0.4, 0.6]. Здесь симметричная относительно х = 0.5 расчетная область не сокращена вдвое по причине желания проверить симметрию вычислительного алгоритма (прогоночные коэффициенты вычисляются слева направо, сеточные функции — справа налево).
5. Результаты расчетов
На рис. 1, а представлено пространственное распределение растепленности с учетом адиабатического расширения газа (13) для моментов времени 1, 10, 50 с. На рис. 1, б приведено также распределение растепленности для случая (12), когда адиабатическое расширение газа отсутствует. В первом случае дополнительного растепления гид-ратизированной среды на границах области действия источников (х Е [0.4, 0.6]) не происходит. Другие параметры (Бт, Р, Т) в обоих случаях ведут себя подобным образом, их пространственные распределения на те же моменты времени представлены на рис. 2.
В случае отсутствия адиабатического расширения газа (т. е. при выборе источников в среде коллектора в виде (12)) проведены дополнительные анализ и тестирование механизмов возникновения пиков растепленности на границе области действия источников (х Е [0.4, 0.6]). Их возникновение связано с тем, что в силу задания отбора в виде (12) с а = 0 вода, образующаяся при разложении гидрата, несжимаема, не отбирается и накапливается на границах зоны стока, что и приводит к соответствующим пикам растепленности при продолжающемся отборе газа на границах области действия зоны стока. Результаты расчетов с разными значениями фазовой проницаемости показали, что нелинейность относительной фазовой проницаемости по воде и газу, а также абсолютной проницаемости не влияет на возникновение пиков. Если отключить отбор после первой секунды, то пики не появляются.
Для оценки практической точности вычислений проводилась серия расчетов с измельчением сетки по времени и пространству от "грубой" сетки с шагами (т\ = 0.01 с и = 0.01 м) до сетки в 16 раз более мелкой (т0 = 0.000625 си = 0.000625 м). Результаты расчетов на самой мелкой сетке брались как эталонные (предполагалось, что они близки к точному решению) и использовались для оценки реальной погрешности и порядка скорости сходимости к "точному" решению. Погрешность оценивалась в норме "С" как максимум модуля разницы сеточных функций, отвечающих решениям на грубой и мелкой сетках. Скорость сходимости оценивалась через логарифм отношения погрешностей, полученных на двух грубых сетках с шагами, различающимися в 2 раза.
а
б
0.66
0.62-
0.58-
0.54-
0.5
0
0.2
0.4
0.6
0.8
0
0.2
0.4
0.6
0.8
Рис. 1. Распределение растепленности для моментов времени 1, 10, 50 с (кривые 1-3 соответственно); а = 0, (3 = 0.00001 при х Е [0.4, 0.6] с учетом (а) и без учета (б) адиабатического расширения газа
Sw 0.680.660.640.620.6
0
0.2
0.4
0.6
0.8
х, м
Рис. 2. Распределение водонасыщенности (а), давления (б) и температуры (в) для моментов времени 1, 10, 50 с (кривые 1-3 соответственно); а = 0, /3 = 0.00001 при х Е [0.4, 0.6] с учетом адиабатического расширения газа
Таким образом, для моментов времени t = 1 и 10 с проводились следующие расчеты: "эталонный" расчет: т0 = 0.000625 с, h0 = 0.000625 м; расчет при п = 1: т\ = 0.01 с, hi = 0.01 м; расчет при п = 2: т2 = 0.005 с, h2 = 0.005 м.
В таблице приводятся полученные результаты погрешностей аппроксимации и порядка скорости сходимости ( Р0 = 10 МПа).
Расчеты проводились на персональном компьютере с тактовой частотой и процессором: 2.7 GHz Intel Core i5. Типичный расчет на сетке из 100 ячеек занимал несколько минут.
а
в
Результаты расчетов
Погрешность аппроксимации Порядок скорости
Параметр при шаге сходимости
п = 1 п = 2 при п = 1, 2
t =1 c
Sw 0.0000011 0.0000005 1.14
Su 0.0000016 0.0000007 1.16
P/Po 0.120839 0.064643 0.90
t = 10 c
Sw 0.020223 0.013281 0.61
Su 0.030380 0.015489 0.97
P/Po 0.062255 0.031878 0.97
Заключение
Сформулирована общая двублочная модель диссоциации гидратов в пористой среде для моделирования теплообмена трех фаз — свободных воды, газа и газогидратов — применительно к полномасштабному предсказательному моделированию на вычислительной технике флюидодинамических и диссипативных процессов в пористых средах, содержащих газовые гидраты в виде твердых включений. Построена разностная схема, позволяющая численно решать систему уравнений фильтрации жидкостей и газов в пористой среде с учетом диссоциации газовых гидратов, неявная по давлению и явная по водонасыщенности и растепленности.
Проведены вычислительные эксперименты, в том числе с учетом адиабатического расширения газа в коллекторном пространстве, результаты которых демонстрируют применимость разработанных методов для расчета реальных задач, связанных с залежами газогидратов, а также с исследованиями сложных процессов динамики водо-и гидратонасыщенности в пласте.
Благодарности. Работа выполнена при частичной поддержке РФФИ (проекты № 16-29-15081-офи_м, 18-07-00841-a) и частично поддержана Проектом повышения конкурентоспособности МИФИ по договору с Министерством образования и науки РФ № 02.A03.21.0005.
Список литературы / References
[1] Истомин В.А., Якушев В.С. Газовые гидраты в природных условиях. М.: Недра, 1992. 236 с.
Istomin, V.A., Yakushev, V.S. Gas hydrates in natural conditions. Moscow: Nedra, 1992. 236 p. (In Russ.)
[2] Гинсбург Г.Д., Соловьев В.А. Субмаринные газовые гидраты. СПб.: ВНИИ океанологии, 1994. 200 с.
Ginsburg, G.D., Solovev, V.A. Submarine gas hydrates. SPb.: VNII Okeanologiya, 1994. 200 p. (In Russ.)
[3] Englezos, P. Clathrate hydrates // J. of Industrial and Eng. Chemistry. 1993. Vol. 32. P. 1251-1274.
[4] Гудзенко В.Т., Вареничев А.А., Громова М.П. Газогидраты. Информационно-аналитический обзор // Геология, геофизика и разработка нефтяных и газовых месторождений. 2016. № 5. С. 39-68.
Gudzenko, V.T., Varenichev, A.A., Gromova, M.P. Gas hydrates. Informational and analytical review // Geologiya, Geofizika i Razrabotka Neftyanykh i Gazovykh Mestorozhdeniy. 2016. No. 5. P. 39-68. (In Russ.)
[5] Басниев К.С., Кочина И.Н., Максимов В.М. Подземная гидромеханика. М.: Недра, 1993. 416 с.
Basniev, K.S., Kochina, I.N., Maksimov, V.M. Underground hydromechanics. Moscow: Nedra, 1993. 416 p. (In Russ.)
[6] Черский Н.В., Бондарев Е.А. О тепловом методе разработки газогидратных залежей // Докл. АН СССР. 1972. Т. 203, № 3. C. 550-552.
Cherskii, N.V., Bondarev, E.A. The thermal method for the development of gas hydrate deposits // Dokl. Akad. Nauk SSSR. 1972. Vol. 203, No. 3. P. 550-552. (In Russ.)
[7] Веригин Н.Н., Хабибулин И.Л., Халиков Г.А. Линейная задача о разложении гидратов газа в пористой среде // Изв. АН СССР: Механика жидкости и газа. 1980. № 1. C. 174-177.
Verigin, N.N., Khabibullin, I.L., Khalikov, G.A. The linear problem of the dissociation of the hydrates of a gas in a porous medium // Fluid Dynamics. 1980. Vol. 15. P. 144-147.
[8] Verigin, N.N., Khabibullin, I.L., Khalikov, G.A. Axisymmetric problem of heat and mass transfer in saturated porous medium // J. of Engineering Physics. 1980. Vol. 38, No. 5. P. 581-585.
[9] Бондарев Е.А., Максимов А.М., Цыпкин Г.Г. К математическому моделированию диссоциаций газовых гидратов // Докл. АН СССР. 1989. Т. 308, № 3. C. 575-578. Bondarev, E.A., Maksimov, A.M., Tsypkin, G.G., Mathematical modeling of gas hydrate dissociation // Dokl. Akad. Nauk SSSR. 1989. Vol. 308, No. 3. P. 575-578. (In Russ.)
[10] Нигматулин Р.И., Ш^агапов В.ШШ., Сыртланов В.Р. Автомодельная задача о разложении газогидратов в пористой среде при депрессии и нагреве // ПМТФ. 1998. Т. 39, № 3. C. 111-118.
Nigmatulin, R.I., Shagapov, V.Sh., Syrtlanov, V.R. Self-similar problem of decomposition of gas hydrates in a porous medium upon depression and heating // J. of Appl. Mechanics and Techn. Physics. 1998. Vol. 39, No. 3. P. 421-427.
[11] Цыпкин Г.Г. Течения с фазовыми переходами в пористых средах. М.: Физматлит, 2009. 232 c.
Tsypkin, G.G. Flows with phase transitions in porous media. Moscow: Fizmatlit, 2009. 232 p. (In Russ.)
[12] Джафаров Д.С. Математическое моделирование диссоциации газогидратов в приложении к интерпретации исследований скважин газогидратных месторождений на нестационарных режимах фильтрации: Дис. ... канд. техн. наук. М.: РГУНГ, 2015. 120 c. Dzhafarov, D.S. Mathematical modeling of dissociation of gas hydrates related to the interpretation of the wells of gas hydrate deposits in non-stationary filtration regimes: Diss. ... kand. tekhn. nauk. Moscow: RGUNG, 2015. 120 p. (In Russ.)
[13] Vasil'eva, Z.A., Efimov, S.I., Yakushev, V.S. Prediction of thermal interaction between oil/gas wells and intra-permafrost metastable gas hydrates // Earth's Cryosphere. 2016. Vol. 20, No. 1. P. 60-63.
[14] Khasanov, M.K., Musakaev, N.G. Gas hydrate formation in porous ice rich methane reservoirs upon injection of carbon dioxide: forward modeling // Earth's Cryosphere. 2016. Vol. 20, No. 3. P. 59-65.
[15] ШШагапов В.ШШ., Чиглинцева А.С., Русинов А.А. Математическое моделирование процесса образования гидрата в пласте насыщенного снегом при нагнетании холодного газа // Вычисл. механика сплошных сред. 2016. Т. 9, № 2. C. 173-181.
Shagapov, V.Sh., Chiglintseva, A.S., Rusinov, A.A. Mathematical modeling of hydrate formation in a reservoir saturated with snow by cold gas injection // Comput. Continuum Mechanics. 2016. Vol. 9, No. 2. P. 173-181. (In Russ.)
[16] Ahmadi, G., Ji, C., Smith, D.H. A simple model for natural gas production from hydrate decomposition // Annals of the New York Acad. Sci. 2000. Vol. 912. P. 420-427.
[17] Kamath, V.A., Godbole, S.P. An analytic model for analyzing the effects of dissociation of hydrates on the thermal recovery of heavy oils // SPE RE. 1988. Vol. 3, No. 2. Paper SPE-14224-PA.
[18] Kamath, V.A., Mutalik, P.N., Sira, J.H. Experimental study of brine injection and depressurization methods for dissociation of gas hydrates // SPE Formation Evaluation. 1991. Vol. 6, No. 4. P. 477-484.
[19] Kim, H.C., Bishnoi, P.R., Heidemann, R.A., Rizvi, S.S.H. Kinetics of methane hydrate decomposition // Chemical Eng. Sci. 1987. Vol. 42, No. 7. P. 1645-1653.
[20] Yousif, M.H., Abass, H.H., Selim, M.S., Sloan, E.D. Experimental and theoretical investigation of methane-gas-hydrate dissociation in porous media // SPE Reservoir Eng. 1991. Vol. 6, No. 1. P. 69-76.
[21] Aziz, K., Settari, A. Petroleum reservoir simulation. London; New York: Elsevier Appl. Sci. Publ., 1979. 476 p.
[22] Goel, N., Wiggins, M., Shah, S. Analytical modeling of gas recovery from in Situ hydrates dissociation // J. of Petroleum Sci. and Eng. 2001. Vol. 29, No. 2. P. 115-127.
[23] Khataniar, S., Kamath, V.A., Omenihu, S.D., Patil, S.L., Dandekar, A.Y. Modeling and economic analysis of gas production from hydrates by depressurization method // Canadian J. of Chemical Eng. 2002. Vol. 80, No. 1. P. 135-143.
[24] Jeannin, L., Bayi, A., Renard, G., Bonnefoy, O., Herri, J.M. Formation and dissociation of methane hydrates in sediments. Pt II: Numerical modeling // Proc. 4th Intern. Conf. on Gas Hydrates, Yokohama, Japan, May 2002. Vol. 2. P. 802-806.
[25] Wilder, J.W., Moridis, G.J., Wilson, S.J. et al. An international effort to compare gas hydrate reservoir simulators // Proc. of 6th Intern. Conf. on Gas Hydrates (ICGH 2008), Vancouver, Canada, July 2008. 12 p.
[26] Pinero, E., Hensen, C., Haeckel, M. et al. 3-D numerical modelling of methane hydrate accumulations using PetroMod // Marine and Petroleum Geology. 2016. Vol. 71. P. 288-295.
[27] Qorbani, K., Kvamme, B. Non-equilibrium simulation of CH4 production from gas hydrate reservoirs through the depressurization method //J. of Natural Gas Sci. and Eng. 2016. Vol. 35. P. 1544-1554.
[28] Vasil'ev, V.I., Popov, V.V., Tsypkin, G.G. Numerical investigation of the decomposition of gas hydrates coexisting with gas in natural reservoirs // Fluid Dynamics. 2006. Vol. 41, No. 4. P. 599-605.
[29] Bondarev, E.A., Rozhin, I.I., Popov, V.V., Argunova, K.K. Assessment of possibility of natural gas hydrates underground storage in permafrost regions // Earth's Cryosphere. 2015. Vol. 19, No. 4. P. 58-67.
[30] Повещенко О.Ю., Гасилова И.В., Галигузова И.И. и др. Об одной модели флюи-додинамики в пористой среде, содержащей газогидраты // Матем. моделирование. 2013. Т. 25, № 10. C. 32-42.
Poveshchenko, O.Iu., Gasilova, I.V., Galiguzova, I.I. et al. One model of fluid dynamics in a porous medium containing gas hydrates // Math. Models and Comput. Simulations. 2013. Vol. 25, No. 10. P. 32-42. (In Russ.)
[31] Казакевич Г.И., Клочкова Л.В., Повещенко Ю.А., Тишкин В.Ф. Математическое исследование системы уравнений газогидратных процессов в пористой среде // Журн. Средневолжского матем. общества. 2011. Т. 13, № 1. C. 7-11.
Kazakevich, G.I., Klochkova, L.V., Poveshchenko, Yu.A., Tishkin, V.F.
Mathematical investigation of a system of equations of gas hydrate processes in a porous medium // Zhurn. Srednevolzhskogo Matem. Obshchestva. 2011. Vol. 13, No. 1. P. 7-11. (In Russ.)
[32] Гасилов В.А., Гасилова И.В., Клочкова Л.В. и др. Разностные схемы на основе метода опорных операторов для задач динамики флюидов в коллекторе, содержащих газогидраты // Журн. вычисл. математики и матем. физики. 2015. Т. 55, № 8. C. 1341-1355.
Gasilov, V.A., Gasilova, I.V., Klochkova, L.V. et al. Difference schemes based on the support operator method for fluids dynamics problems in a collector containing gas hydrates // Comput. Math. and Math. Phys. 2015. Vol. 55, No. 8. P. 1310-1328.
Поступила в 'редакцию 29 марта 2018 г.
The model of fluid dynamics in a hydrate-containing porous medium with splitting into physical processes
Rahimly, Parvin I.1, Sharova, Yulia S.2'*, Rahimly, Orkhan R.1, Podryga, Viktoriia O.2'3, Gasilova, Irina V.2, Popov, Sergey B.2, Poveshchenko, Yuriy A.2'4
1Moscow Institute of Physics and Technology, Dolgoprudny, Moscow Region, 141701, Russia
2Keldysh Institute of Applied Mathematics of RAS, Moscow, 125047, Russia
3 Moscow Automobile and Road Construction State Technical University, Moscow, 123182, Russia
4National Research Nuclear University MEPhI, Moscow, 115409, Russia * Corresponding author: Sharova, Yulia S., e-mail: [email protected]
Goal. It is necessary to understand the problems of the phase transformations of gas hydrate inclusions in the porous media since it can help in the development of the technologies of gas hydrate production.
Method. The paper deals with the thermodynamically equilibrium model of the splitting into physical processes of a two-component three-phase filtration fluid dynamics with gas hydrate inclusions, for which a family of two-layer completely conservative difference schemes based on support operators method with space-time temporal scales is constructed. The initial problem, formulated in the form of the first principles (conservation of H2O, CH4 mass and total energy of the medium), with a common matrix of the system relative to the thaw, moisture saturation, pressure and temperature, has mixed hyperbolic and multiscale parabolic properties. The direct unsplitted utilization of such system for the purpose of determining the dynamics of variables and constructing an implicit difference scheme, which is required for calculations of filtration parabolic processes with large steps in time, is difficult. This model is based on splitting into physical processes.
Result. The numerical model allows discretization of the problem in an one-dimensional case and implementing unconditionally stable difference scheme. The obtained results demonstrate the applicability of the proposed model for the solution of typical problems of gas hydrate fluid dynamics, including the studies of the complex dynamics of water and hydrate saturation of the formation in respect of adiabatic expansion of the gas in the collector space.
Value. The application of these schemes is also important in the combination modeling of the three-phase medium with processes in the thawed zone where the
© ICT SB RAS, 2018
hydrate is absent, since the number of unknowns, equations, their types here change, and it is important at the difference level to transfer the balances divergently into different flow areas. The approach outlined in the work is the basis for constructing algorithms for transphase hydrate encapsulation.
Keywords: gas hydrates, filtration, mathematical modeling.
Cite: Rahimly, P.I., Sharova, Yu.S., Rahimly, O.R., Podryga, V.O., Gasilova, I.V., Popov, S.B., Poveshchenko, Yu.A. The model of fluid dynamics in a hydrate-containing porous medium with splitting into physical processes // Computational Technologies. 2018. Vol. 23, No. 6. P. 64-79. (In Russ.) DOI: 10.25743/ICT.2018.23.6.007.
Acknowledgements. This research was partially supported by RFBR (projects No. 16-29-15081-ofi_m, 18-07-00841-a). This work was partially supported by the Academic Excellence Project of the NRNU MEPhI under contract with the Ministry of Education and Science of the RF No. 02.A03.21.0005.
Received 29 March 2018
