Математическая модель и алгоритм решения задачи неизотермической фильтрации газа в пласте с учетом разложения гидрата Текст научной статьи по специальности «Математика»

УДК 536.2, 519.63

DOI: 10.14529/mmph170203

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И АЛГОРИТМ РЕШЕНИЯ ЗАДАЧИ НЕИЗОТЕРМИЧЕСКОЙ ФИЛЬТРАЦИИ ГАЗА В ПЛАСТЕ С УЧЕТОМ РАЗЛОЖЕНИЯ ГИДРАТА

Н.Г. Мусакаевл, С.Л. Бородин1, Д.С. БельскиХ

1 Тюменский филиал Института теоретической и прикладной механики им. С.А. Христиановича СО РАН, г. Тюмень, Российская Федерация

2 Тюменский государственный университет, г. Тюмень, Российская Федерация E-mail: bedeser@yandex.ru

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

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

Введение

В настоящее время углеводороды являются одним из важнейших источников для энергетики и практическим незаменимым сырьем для химической промышленности. В связи с этим ведутся исследования в области поиска и разработки новых, нетрадиционных источников углеводородов, одним из которых особый интерес представляют гидраты природного газа [1, 2]. Запасы углеводородов в гидратах огромны - специалисты сходятся на том, что в скоплениях природных газовых гидратов содержится примерно 20 000 трлн. м3 метана (для примера в традиционных извлекаемых запасах - 250 трлн. м3) [3]. Таким образом, задачи, связанные с добычей метана из газовых гидратов являются чрезвычайно перспективными.

Для решения практических задач требуется теоретическая проработка, основывающаяся на адекватных математических моделях и численных методах их реализации. В настоящей работе предложена математическая модель неизотермической фильтрации газа, а также разработан алгоритм решения уравнений этой модели, который включает в себя авторский метод расчёта гидратонасыщенности. Наибольший интерес представляет именно расчет гидратонасыщенности, так как уравнения неизотермической фильтрации газа без учета фазовых переходов решаются относительно легко, а учет фазового перехода создает значительные вычислительные сложности. Например, одним из авторов проводилось сравнение численных методов решения задачи Стефана о фазовом переходе [4], и ни один из семи рассматривавшихся в этой работе методов не подходит для расчета процесса разложения/образования гидрата в пласте. Поэтому предложен метод, который позволяет решать такую задачу, причем он является универсальным - подходит и для многомерных, и для многофронтовых задач Стефана, а также для задач с протяженной областью фазовых переходов.

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

Построим математическую модель неизотермической фильтрации газа с учетом фазового перехода «гидрат метана « метан + вода». При этом будем рассматривать нагнетание газа (метана) в однородный горизонтальный пористый пласт постоянной толщины. Данная залежь в исходном состоянии заполнена газом (одноименный закачиваемому) и его гидратом, давление p0 и температура T0 которых соответствуют термодинамическим условиям их существования в свободном состоянии. Температура нагнетаемого в пласт газа Te является постоянной величиной и Te > T0.

Примем следующие допущения: гидрат является двухкомпонентной системой с массовой концентрацией газа G; кровля и подошва пласта непроницаемы; температуры газа, жидкости,

гидрата и пористой среды в каждой точке пласта совпадают (однотемпературная модель). Положим, что пористость пласта т постоянна, скелет пористой среды, газовый гидрат и вода несжимаемы и неподвижны [2, 5].

Запишем уравнения сохранения массы для газа, воды и гидрата [2]:

Э

дг

(т8гРг) + (т8гР^г) = ^ = GJ : ^(Ръ ) = = (1 - 6)J, ^ (mShРh ) = -Jh®g - =^,

$ ъ + Sh + Sg = 1 ,

(1) (2)

(3)

(4)

где $ (7 = g, ъ, И) - насыщенность поры 7-ой фазой; р7 (7 = g, ъ, И) - плотность 7-ой фазы; vg - скорость газа, J - интенсивность разложения гидрата; параметры газа, воды и гидрата снабжены индексами g, ъ и к соответственно.

Скорость фильтрации газа удовлетворяет уравнению Дарси

mSgVg = —-(8гаа р - р^), (5)

mg

где - проницаемость по газу; - вязкость газа; р - давление; g - ускорение свободного падения.

Для газа примем следующие уравнения состояния и коэффициента сверхсжимаемости г [6]:

Р = гряКГ , (6)

(

г =

0,418

( г ^ Г

V с

\

+ 0,73

Рс р

+ 0,1—.

(7)

где Я - удельная газовая постоянная; Г - температура; рс, Гс - критические давление и температура газа.

Запишем уравнение энергии в общем виде [7] с учетом фазового перехода:

^ ~ ^ (1 лтЛ т

дГ

рс +рgcgmSgv

grad Г + е grad р--

cg )

PgcgmSghs — = Шу (Г) + трь1^—-

е =-

д(

(8)

1 г (дг =

Pgcg г V дГ )р' - Pgc:

-е

g g

рс = (1 - т)рксЯк + т (Sgрg +ShPhCh),

1 = (1 - т)Лк + т (SgЛg + ^Лъ + ^л- ).

Здесь Ьк - теплота разложения газового гидрата; е7 - коэффициент Джоуля-Томпсона, п - коэффициент адиабатического охлаждения; р-к - плотность скелета пористой среды; с7 и Х7 (7 = -к, g, ъ, к) - удельная теплоёмкость и коэффициент теплопроводности скелета, газа, воды и гидрата, соответственно.

Рассмотрим задачу в одномерной плоскопараллельной постановке. Тогда после преобразований уравнений (1)-(6) и (8) получим следующую систему:

др р дГ р дг р д7 _ Г д7 г д^ ^ д^

д^?^ гЯГ д

дТ = рмс1_ д7 рс

(

+---

Sg дх

2

рёкг

тт

\

др

дх

- .р. д^

^др дГ + е ^( др? + mSgrls др Ц дх дх V дх) а д1

рк

+Xд(лэг)+тр-ьк ^к

рс дх V дх) рс д7

^ = ^0 + —(1 - 6)(Sh0 - Sh ). ръ

(9)

(10) (11)

с

Зададим начальные и граничные условия для задачи закачки теплого газа в пласт, изначально содержащий газ и его гидрат:

г = 0, Xw < х < Ь: р = Ро, т = T0, 8к = Sh0, ^ = 0, ^ = 1 - ^ (12)

Р = Ре, т = те;

t > 0:

х = х„

х = L: дР = 0, ^ = 0. Эх Эх

(13)

Здесь хн, - радиус скважины, Ь - длина пласта, р0, Т0 - начальные давление и температура в пласте, Sh0 - начальная гидратонасыщенность, ре - давление на забое скважины, Те - температура закачиваемого газа.

Уравнения (9)-(11), уравнение Латонова-Гуревича для сверхсжимаемости (7) и уравнение (4) для расчета газонасыщенности образуют систему, которая позволяет рассчитывать распределения давления, температуры, насыщенностей и других параметров в пласте при неизотермической фильтрации газа, которая сопровождается разложением или образованием гидрата. В данной системе 5 уравнений и 6 неизвестных, поэтому для ее замыкания используется описанный далее численный метод решения, который заключается в расчете гидратонасыщенности, исходя из следующих лимитирующих условий:

• при разложении гидрата давление повышается, пока не достигнет равновесного давления или не разложится весь гидрат;

• при разложении гидрата температура понижается, пока не достигнет равновесной температуры или не разложится весь гидрат;

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

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

Данные условия служат для уточнения итерационных параметров, которые подставляются в уравнения (4), (7), (9)—(11) до тех пор, пока не будет достигнута приемлемая погрешность. Алгоритм решения подробно описан далее.

Алгоритм решения

В работе при численной реализации математической модели были использованы неявная разностная схема, метод прогонки, метод простых итераций и новый метод для расчета гидратонасыщенности. Дискретные аналоги уравнений (9) и (10) могут быть представлены в следующем виде:

/ \к+1 / \к

(^), -(^) ^ • я• тк+1 1

„k+1 „k „к+1 rjik+1 rpk к+1 k+1 k „..

p, - p, _Pi T -T , Pi Z - z p,

Dt

x

rA

mm

T

k+1

k+1

Dt

- + -

k+1

Dt

k+1 (Sg),

Dt

- + -

k+1 (Sg),

k+1

„k+1 „k+1 ( r k ~\k " jk+1 „k+1 Pgkg

J+1 2

T

Tk

p,+1 - p,

Ах

pgcg

g g mms

pi

■p, -1

2

Ах

k+1

g k+1

Ах

Ph _G (s, )r1 -(s, )k.

x

p,

k+1 k (sg), (Pg),

k +1

Dt

k+1 rrk ( P c \k+1 ( k \k+1 ( - k+1 „k+1 ^

Dt

pc

+

PgC.

\

Ji k+1

V mg J

p,+1 - pi

Ах

jik+1 jik+1 1 i+1 -1 ■

Ах

l— + (e,):

k+1

( pi+11 - p/ Ах

k+1

+

k+1 (h )k+1 pk+1 - p, + mphLh (S, )k+1 -(S,)k +

g g PC

+

m (Sg ), (hs ), 1 1

Dt

(pc )k

k+1

Dt

(pc )k+1 Ах

rpk+1 rpk+1 rpk+1 rpk+1

ik+1 ' i+1 ' i о k+1 ' i ' i-1

1 1--1 1 -

i+- Ах i— Ах

22

(14)

(15)

Здесь верхние индексы относятся к шагу по времени, а нижние - по пространству. Для решения системы с помощью метода прогонки необходимо уравнения (14) и (15) привести к специальному виду

Ар • рй1 + Ср • рк+1 + Вр • рк+ = р, Аг • Г-+1 + Сг • Гк+1 + Вг • Г++1 = р , где коэффициенты А, В/, С/ (/ = р, Г) выражаются следующим образом:

ар =

At

гк+1 • Я • Гк+1 рек.

/ Лк+1 к

(Ах )2 ()к+

тт

вр =

А

*к+1 • Я • Гк+1 Г ргк,

\к+1

грк

С =——

р гк+1

к+1

+

к+1 п т^к+1

г • Я • Г

(^ )к _а__

& )к+1 (Ах)2 &)к+1

(* )2 ( рА

к+1

g g

( П Ьг \к +1 ( п Ьг \к +1

тт

+

2

тт тт

/7+-

),-1

рнО

к+1 к+1 (sg), (р, )

(Sh )к+1 -(Sh )к

р =-рк •

* р Г1 5

Аг —

А!

Лк+1

7--

2

(ах) (рс)г

а!

к+1

Вг =

А!

(Ах )2

лЛ+1

7+2 к+1

+

( р с Лк+1 Г к Лк+1 р,с, к■

(рс )к

рс

V т)

( рк+11 - рк+1)

Сг =-1 -

р = -Гк - I ^

(Ах )2

лк+1 Г к Л

к+1

g g рс

/ 7

( рк+11 - рк+1 )

+ -

(рс )к

к+1

Л

Лк +1 +Лк +1

7-- 7+—

V 2 2 )

чк+1

рс

к

(е)к

к +1

( Ь \

к+1

А!

т) (Ах)

(sh )к+1 -(sh )к

(рс ),

( /£1 - рк+1 )2 + т ( Sg )к+1 (1- )к+1 ( рк+1 - рк )

к +1

А!

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

а,

7+1

-В,

Р - А/ я!

А/а + С,

р/ =- } }

,Р 7 +1

А/•а} + С,

где/ = р, Г при этом:

Гк+1 = а • гк+1 + Яг рк+1 = ар • рк+1 + Яр

Г7 =а7+1 Г 7+1 + Я7+1, р7 =а7+1 р7+1 + Я7+1.

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

к = 0, Яг = Г, к+1 = -Я „к+1 = -Яр 1

К = 0, рр = ре.

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

* а -1'

РN

а

к+1

г =

Г

к+1

Л рк +7рс

0,4 • к-7— + 0,73 Г

N

к+1

+ 0,1 •

р,

(16)

Порядок численного решения следующий:

1. Методом прогонки находятся значения давления в узлах пространственной сетки на новом временном шаге, используя уравнение (14). При этом используются итерационные значения ко-

2

2

г

2

1

g

g

1

с

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

2. Рассчитывается значение коэффициента сверхсжимаемости газа из соотношения (16), используя определенные в 1-ом пункте значения давления и итерационные температуры.

3. Методом прогонки вычисляются значения температуры в узлах пространственной сетки на новом временном шаге из уравнения (15). При этом используются значения давления, найденные в пункте 1, коэффициента сверхсжимаемости - в пункте 2 и итерационные значения гидрато-, водо- и газонасыщенности. Параметры, которые зависят от температуры (например, плотность газа), определяются на основе итерационного значения температуры.

4. Для расчета гидратонасыщенности, во-первых, вычисляются значения равновесной температуры и давления фазового перехода для уже найденных в пунктах 1-3 величин давления, коэффициента сверхсжимаемости и температуры, используя корреляцию Хафермана-Миллера [8]

T = Т

B

f

A - ln

f к+1 к+1 ^ Р, • z

106

> Ph = exp

A -

B

T

к+1

106

к+1

где размерность давления - МПа; размерность температуры - К; Ть, рь - равновесные температура и давление; Аь В\ - эмпирические коэффициенты, имеющие различные значения для двух диапазонов температур, в зависимости от того, на какие фазы разлагается (из каких фаз образуется) гидрат (газ и лед, либо газ и вода):

250 < Т < 272,95: А1 = 8,486; В1 = 2 087,40;

272,95 < Т < 287: А1 = 25,659; В1 = 6 774,09.

Если термодинамические условия допускают разложение гидрата, т. е. (Т] )п™ > Ть и гидрат есть ((Sh)/ ег > 0), то рассчитываем максимальное количество гидрата, которое может разложиться:

ASh = min

DST =(Tk+1 - Th , DShp =( ph

V ' mPhLh V

к+1

/ \iter / \k к+1\ ( g )i \PS )i jc, \iter

'p, ) ^73+1 , (Sh)i

GPhPk

тогда гидратонасыщенность на новом шаге

(Sh )*+1 =( Sh Г .

Если термодинамические условия допускают образование гидрата, т. е. (Т )п™ < Ть и есть гидратообразователи (^и,)/'ег > 0 и ^)/'ег > 0), то рассчитываем максимальное количество гидрата, которое может образоваться:

\ ыет / \ к+1

~срнрк

ASh = min

AST =(Th - Tk+1 , DSp =( pk+1 - Ph)(Sg} (Pg+1)' , AS? = ( Г P

1 - G Ph

тогда гидратонасыщенность на новом шаге:

(Sh )к+1 =(Sh + .

5. Рассчитываются значения водонасыщенности из дискретного аналога уравнения (11):

)к + Р

Ръ

6. Рассчитываются значения газонасыщенности, используя уравнение (4):

(^ )к+1=1 -(Sw )к+1 -(Sh )к+1.

7. Определяется относительная погрешность итерации:

( St )к+1 =( ST )к +P (1 - G )(( Sh )к -( Sh )к+1 ) .

8 = тах

V/

Новое итерационное приближение:

рГ - рк+1 иег к+1 г — 2- грыег грк+1 Т / Т /

Р,к+1 к+1 2 Тк+1 Т

/ \ пег / \ к+1 ( 5, ), -( 5, ),

к+1 ( ),

. „^ег „к+1. иег к+1. т<иег грк+1.

V/: р, = р, , г, = гг , Т,/ = Т,, ,

(^ Г =( )к+1; (Г =( ^ )к+1; (^ )Г =( 5^ )к+1.

% ¡1 \ , /

Если погрешность итерации 3 > е (где е - малая величина задающая точность), то итерации продолжаются, т. е. переходим к пункту 1, в противном случае - к пункту 8.

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

Значения теплоёмкости газа, его теплопроводности и вязкости в зависимости от температуры и давления можно найти линейной интерполяцией табличных данных из справочников [9, 10]. Текущая плотность газа будет известна через уравнение состояния газа (6).

Результаты расчетов давления, температуры и гидратонасыщенности представлены на рисунке.

Распределение по координате х давления (а), температуры (б) и гидратонасыщенности (в) в различные моменты времени

Заключение

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

Литература

1. Макогон, Ю.Ф. Гидраты природных газов / Ю.Ф. Макогон - М.: Недра, 1974. - 208 с.

2. Шагапов, В.Ш. Динамика образования и разложения гидратов в системах добычи, транспортировки и хранения газа / В.Ш. Шагапов, Н.Г. Мусакаев. - М.: Наука, 2016. - 238 с.

3. Коллетт, Т.С. Растущий интерес к газовым гидратам / Т.С. Коллетт, Р. Льюис, Т. Учида // Нефтегазовое обозрение. - Осень, 2001. - С. 38-53.

4. Бородин, Л.С. Численные методы решения задачи Стефана / Л.С. Бородин // Вестник Тюменского государственного университета. Физико-математическое моделирование. Нефть, газ, энергетика. - 2015. - Т. 1, № 3(3). - С. 164-175.

5. Мусакаев, Н.Г. Математическое моделирование процесса добычи газа из газогидратной залежи с учетом образования льда / Н.Г. Мусакаев, М.К. Хасанов // Вестник Тюменского государственного университета. Физико-математическое моделирование. Нефть, газ, энергетика. -2014. - № 7. - С.43-50.

6. Латонов, В.В. Расчёт коэффициента сжимаемости природного газа / В.В. Латонов, Г.Р. Гу-ревич // Газовая промышленность. - 1969. - № 2. - С. 7-9.

7. Басниев, К.С. Подземная гидромеханика / К.С. Басниев, И.Н. Кочина, В.М. Максимов. -М.: Недра, 1993. - 416 с.

8. Истомин, В.А. Газовые гидраты в природных условиях / В.А. Истомин, В.С. Якушев. - М.: Недра, 1992. - 236 с.

9. Варгафтик, Н.Б. Справочник по теплофизическим свойствам газов и жидкостей / Н.Б. Вар-гафтик. - 2-е изд., перераб. и доп. - М.: Наука, 1972. - 720 с.

10. Варгафтик, Н.Б. Справочник по теплопроводности жидкостей и газов / Н.Б. Варгафтик, Л.П. Филиппов, А.А. Тарзиманов, Е.Е. Тоцкий. - М.: Энергоатомиздат, 1990. - 352 с.

Поступила в редакцию 27 января 2017 г.

Bulletin of the South Ural State University Series "Mathematics. Mechanics. Physics" _2017, vol. 9, no. 2, pp. 22-29

DOI: 10.14529/mmph170203

MATHEMATICAL MODEL AND ALGORITHM FOR SOLVING THE PROBLEM OF NON-ISOTHERMAL GAS FILTRATION IN RESERVOIR IN CASE OF HYDRATE DECOMPOSITION

11 2 N.G. Musakaev1, S.L. Borodin1, D.S. Belskikh2

Tyumen Branch of Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Tyumen, Russian Federation

2 Tyumen State University, Tyumen, Russian Federation E-mail: bedeser@yandex.ru

The paper formulates a problem of injection into porous bed, filled up in the initial condition with hydrate and gas, warm (with the temperature higher than the initial temperature of the bed) gas. A mathematical model of non-isothermal gas filtration in case of gas hydrate dissociation is developed to solve this problem. The article presents a solution algorithm, where an implicit difference scheme, a sweep method and a method of simple integration are applied. The method for calculating hydrate saturation from several limiting conditions is suggested. It can be used for solution of other phase-change problems, also for multidimensional Stefan problems, as well as problems with an extended phase transition zone. After that the problem is considered in one-dimensional plane-parallel formulation with regard to required initial and boundary conditions for finding a computational solution of a set of equations describing this model. At the end, the paper presents the problem calculation results using the suggested method, on the basis of which the distribution of parameter values for some time intervals are shown. In the performed calculations the reservoir in the initial condition is filled up with methane and its hydrate.

Keywords: non-isothermal filtration of gas; gas hydrate; numerical method; phase transition.

References

1. Makogon Yu.F. Gidraty prirodnykh gazov (Natural gas hydrates). Moscow, Nedra, 1974, 208 p. (in Russ.).

2. Shagapov V.Sh., Musakaev N.G. Dinamika obrazovaniya i razlozheniya gidratov v sistemakh dobychi, transportirovki i khraneniya gaza (Dynamics of hydrate formation and decomposition in the systems for gas production, transportation and storing). Moscow, Nauka Publ., 2016, 238 p. (in Russ.).

3. Kollett T.S., L'yuis R., Uchida T. Neftegazovoe obozrenie, Autumn, 2001, pp. 38-53. (in Russ.).

4. Borodin L.S. Tyumen State University Herald. Physical and Mathematical Modeling. Oil, Gas, Energy, 2015, Vol. 1, no. 3(3), pp. 164-175. (in Russ.).

5. Musakaev N.G., Khasanov M.K. Tyumen State University Herald. Physical and Mathematical Modeling. Oil, Gas, Energy, 2014, no. 7, pp. 43-50. (in Russ.).

6. Latonov V.V., Gurevich G.R. Gazovaya promyshlennost', 1969, no. 2, pp. 7-9. (in Russ.).

7. Basniev K.S., Kochina I.N., Maksimov V.M. Podzemnaya gidromekhanika (Subsurface hydromechanics), Moscow, Nedra Publ., 1993, 416 pp. (in Russ.).

8. Istomin V.A., Yakushev V.S. Gazovye gidraty v prirodnykh usloviyakh (Gas hydrates in natural environment), Moscow, Nedra Publ., 1992, 236 p. (in Russ.).

9. Vargaftik N.B. Spravochnikpo teplofizicheskim svoystvam gazov i zhidkostey (Reference book on heat-transfer properties of gas and liquid). Moscow, Nauka Publ., 1972, 720 p. (in Russ.).

10. Vargaftik N.B., Filippov L.P., Tarzimanov A.A., Totskiy E.E. Spravochnik po teploprovodnosti zhidkostey i gazov (Reference book on heat conductivity of liquid and gas). Moscow, Energoatomizdat Publ., 1990, 352 p. (in Russ.).

Received January 27, 2017