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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ

УДК 519.6:532

О. С. БОРЩУК

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

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

ВВЕДЕНИЕ

Задача полномасштабного гидродинамического моделирования процессов нефтегазодобычи сопряжена с ресурсоемкими вычислениями, кроме того, размеры моделей постоянно увеличиваются. Так, еще 5-6 лет назад большинство гидродинамических моделей содержало не более 100 тысяч ячеек, теперь же средний размер модели достигает 1 млн ячеек и более. С учетом роста размерности моделей особую роль играет эффективность методов решения разреженных систем линейных алгебраических уравнений большой

размерности, тем более, что процент времени решения линейной системы в общем времени расчета составляет от 80% до 90%.

В статье предложена модификация предобуславливателя ОРЯ, который совместно с каким либо методом подпространства Крылова, является одним из самых быстрых при решении линейных систем в пакетах гидродинамического моделирования залежей углеводородов. Тестирование нового метода проводилось на 2-х и 3-х фазных моделях и показало преимущество в скорости в среднем на 20%, в количестве требуемой оперативной памяти на 35-40%.

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим систему уравнений, описывающих двухфазную фильтрацию вязкой сжимаемой жидкости в пористой среде [1] (все

изложенное ниже может быть продолжено на случай трехфазной фильтрации):

ам„

дґ

м

дґ

-У(КХ № ) =

-У{КХ0УФ0 ) = д0,

(1)

(2)

где М, = фрЛ, Мо = ФРо (1- Sw), I р =р р -

Ц р

мобильность фазы р={о, w}, Фо = р -Роё^Е ,

Ф * = Р + Рсош - Р , * - врем^ ф = ф(р) -

пористость, Sw - насыщенность воды, р, = р,(р, Sw), Ро = р0(р) - плотности воды и нефти,

к =

УУ

к.

- диагональный тензор

абсолютной проницаемости, к*, = к,,^,), кго = = кго^„) - относительные фазовые проницаемости воды и нефти, Ц, = Ц,(р, Sw), Ц0 = Цо(р) -вязкости воды и нефти, р - давление, рС0№= = рС0ш (Я„) - капиллярное давление между

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

Отметим, что уравнения (1, 2) представляют собой уравнения сохранения массы для воды и нефти.

Дискретизируя уравнения (1, 2) по

пространству с помощью метода конечных

объемов, получим систему нелинейных

уравнений:

V (к) (м*) = £ Р, + еН, (3)

]'еМ1 (Н)

уг (к)-(Мо) = £ Уо] + ~@о, V/ е Н, (4)

]еЫ‘ (Н)

где О — множество ячеек дискретной сетки, У - объем г-ой ячейки, N(0.) - множество соседних ячеек для ячейки г, М - шаг по времени, Q - источник или сток, а - означает, что значение переменной (функции) а берется с нового временного слоя,

Р, = Т1,]к, (ф* -Ф*) - поток водной фазы через границу между ячейками г и ], рч = т1’}хо (фо -фго) - поток нефтяной фазы через границу между ячейками г и , Тг, -проводимость между ячейками г и ] (зависит от тензора К ), ё - ускорение свободного падения. Отметим, что мобильности к*] и к о вычисляется против потока.

Нелинейная система (3, 4) состоит из 2 X п уравнений (п - число ячеек сетки) и 2 X п неизвестных (давление и насыщенность воды в каждой ячейке). Для решения систем нелинейных алгебраических уравнений наиболее часто используется метод Ньютона. Возникающая при этом несимметричная линейная система алгебраических уравнений (матрица Якоби) представляет собой разреженную матрицу состоящую из блоков:

ЭRL ЭRi

Э^

R

ЭS^

Эp1

R

ЭpJ

i, є W,

(5)

где - уравнения (3), а Яо - уравнение (4).

Отметим, что в строке 1 линейной системы количество ненулевых блоков равно количеству соседей ячейки г плюс диагональный блок.

2. ПРЕДОБУСЛАВЛИВАТЕЛЬ СРЯ

Как было показано ранее, система линейных алгебраических уравнений

Ах = Ь, (6)

получающаяся в методе Ньютона при решении системы (3, 4) представляет собой несимметричную разреженную матрицу с блочной структурой. При этом все блоки имеют размер 2 X 2 .

Наиболее популярными и быстрыми методами решения разреженных систем большой размерности являются методы подпространства Крылова GMRES [4] и BiCGStab [4], в которых используются различные методы предобуславливания. При этом вместо системы (6) решается эквивалентная система

(7)

где M - матрица предобуславливателя, удовлетворяющая следующим условиям:

1. det(M) ф 0;

2. MT1 - легко вычисляема;

3. M ~ A.

Отметим, что скорость сходимости метода напрямую зависит от выбора предобуслав-ливателя M.

Сложность решения рассматриваемых нами систем связана со смешанным характером исходной системы уравнений (1, 2). Так, в [1] показано, что система (1, 2) является

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

дифференциальных уравнений в частных производных соответствующего типа.

Известно, что алгебраический

многосеточный метод (AMG) [З] хорошо

подходит (как предобуславливатель) для решения эллиптических и параболических линейных систем, но при этом неэффективен для гиперболических систем. С другой стороны неполное LU разложение (ILU0) [4] и его модификации хороши при решении

гиперболических систем, но с трудом справляются с эллиптическими и

параболическими системами.

Так как линейная система имеет смешанный характер, был предложен двухступенчатый предобуславливатель Constrained Pressure Residual

(CPR) [2]:

• из исходной линейной системы построить систему на давление при условии, что насыщенность меняется слабо;

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

• полученное решение на давление использовать как начальное приближение для неполного LU разложения (вторая ступень CPR).

Таким образом, CPR позволяет эффективно решать рассматриваемую систему смешанного типа. Результаты численных экспериментов можно найти в [5].

В матричном виде двухступенчатый предобуславливатель может быть записан как

М- = М 2-1 [Е - ЛМ~1 ]+М-1, (8)

где М1 и М2 - матрицы предобуславливателей первой и второй ступени, Е - единичная матрица. В случае предобуславливателя CPR (8) примет вид:

М-РК = М~1и [Е - ЛМр ]+Мр.

Рассмотрим алгоритм метода CPR подробнее. Для построения системы на давление отсортируем столбцы и строки линейной системы (6) таким образом, чтобы вначале шли переменные на давление (индекс р), а потом переменные на насыщенность (индекс 5):

(9)

Предполагая, что насыщенность в ячейках меняется слабо и не зависит от соседних ячеек, исключим из Лрх все недиагональные элементы ВРх=В(ЛРх), то же самое проделаем и с Л^, то есть 06.6=0(Л6.6). С учетом этого приближения система (9) примет вид

Пусть нам надо решить MCPRx = b (или

1 A pp I ps A 1 x p I 1 о p I

I Asp I Ass 1 sx 1 I bs I

" App I ps Dp 1 x p I I p b 1

_ Asp I Dss I xs I I bs I

(lO)

Далее, так как - диагональная матрица, а

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

х^ из второго матричного уравнения системы (10)

Г.-1

xs=Dlbs - DssAspxp

(ll)

и подставив это выражение в первое матричное уравнение системы (10) получим систему на давление:

(Лрр - вр*в- Ар)хр=ър- вр°-ь*’

или

лрхр=гр .

В матричном виде систему на давление можно представить как Лр = СтЫЛС, где

E - DpsD-sl

CT = [E O], N = ps ss , a матрицу

I O E I

предобуславливателя Mp как M-1 = CA-lCTN .

x = M-1pRb), тогда алгоритм будет следующим:

1) rp=bp - DpsD-lbs ( rp=CTNb ) - расчет

правой части для системы на давление;

2) Ap=App - DpsD-sl Asp (Ap=CTNAC) -

построение системы на давление;

3) Ap5xp = rp - решение системы на дав-

ление;

4) 5x = (6xp,0)T (Sx = C6xp ) - расширение вектора Sxp ;

5) b = b - A6x - обновление правой части;

6) MILUx = b - решение системы, где Milu=ILU 0(A) матрица неполного LU разложения для A;

7) x = x + Sx - коррекция решения. Отметим, что в пункте 4 алгоритма не используется выражение переменных xx через переменные xp, это связано с нестабильностью получаемых при этом результатов вблизи источников, стоков или фронта вытеснения.

3. ПРЕДОБУСЛАВЛИВАТЕЛЬ ACPR

При всех преимуществах предобуславливателя CPR (быстрая сходимость, фактически линейный рост времени решения в зависимости от размерности задачи) он обладает одним существенным недостатком. Этим недостатком является повышенное требование к объему оперативной памяти. Так, для решения системы с использованием предобуславливателя CPR требуется в среднем в 1,5-2,0 раза больше оперативной памяти, чем при использовании предобуславливателя ILU0.

В данном разделе предложена модификация предобуславливателя CPR, названная ACPR (Adaptive Constrained Pressure Residual), которая позволяет сократить затраты оперативной памяти, а также уменьшить общее время построения решения.

Основная идея ACPR состоит в разделении всего множества ячеек сетки Q на два непересекающихся подмножества Qe и Q;, где

- множество ячеек, для которых решение на насыщенность может быть получено из решения на давление с использованием формулы (11), а Q = Q/Q е - все оставшиеся ячейки. Так как решение на давление мы получаем на первом шаге предобуславливателя CPR, то для ячеек из множества Qe строить неполное LU разложение не требуется. Таким

i,j

O;

образом, матрица предобуславливателя М^и строится не на основе матрицы А, а на основе матрицы А , которая получается из матрицы А по следующим правилам:

• если г е Ог и у е Ог , то агу = а,

• если г е и у е Оа , то а; .

I ^ а I,]

• если г = ] и г е Ое , то аг у = аг у ;

• если г е Ое и г Ф ], то аг у = 0.

Здесь I - индекс строки, j - индекс столбца, аг у и аг у - блоки матриц А и А

соответственно.

Для разделения множества О на множества

и Ое будем использовать критерий CFL (Courant-Friedrichs-Lewy) [6] (в русскоязычной

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

CFL = max(CFLw , CFLio),

(іг)

где

(

CFLiw = AtApi l'w • max

к^кж

Axi ’ Ay

д.

Az1

CFL0 = AtAp' X0 • max

kL k' hi

\

_______yy___

Ax’ ’ Ay1 ’ Az1

Если выполнено условие CFL1 < 1, то ячейка i является устойчивой.

Будем выбирать множество Qe по

следующему алгоритму:

1) Qe = Q,

2) Qe = Qe/{i: | QW \ + \Ql0\> 0} - убрать все ячейки в которых есть источник или сток,

3) Qe = Qe/{i: CFL > 1} — убрать все ячейки для которых не выполнено условие CFL.

4. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

Численные результаты получены в пакете гидродинамического моделирования трехфазных пластовых систем углеводородов «BOS» нефтяной компании Роснефть.

При тестировании были использованы следующие гидродинамические модели:

1) двухфазная модель SPE10 — 1122000 ячеек, 4 добывающих скважины, 1 нагнетательная;

2) трехфазная модель - около 850000 ячеек, 83 добывающих скважины, 17 нагнетательных.

Полученные результаты представлены в таблице (1)

Т аблица 1

Результаты замеров, где N - число временных шагов, Ы„ - число Ньютоновских итераций, N -число линейных итераций, Т - общее время решения линейных систем, Ут - объем оперативной памяти

Модель 1 2-х фазная Модель г 3-х фазная

Nt 408 1б3

Nn 1235 540

CPR Ni 3879 г371

T (минут) 348 403

Vm (ГБ) 1,25 1,89

ACPR Ni 3886 г399

T (минут) 276 гб8

Vm (ГБ) 0,83 1,08

При

Qe

этом доля ячеек множества ое относительно множества О составила для 2-х фазной модели 3%, а для 3-х фазной 4,5%.

Отметим, что все результаты

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

|Ь - Ах\\

„ „ 2 < 10-4.

ВЫВОДЫ

Предложена модификация предобуслав-ливателя СРЯ, позволяющая сократить время решения разряженных систем линейных алгебраических уравнений, возникающих при гидродинамическом моделировании залежей углеводородов, на 20% и объема требуемой оперативной памяти на 35-40%. Проведено тестирование предложенного алгоритма на 2-х и 3-х фазных задачах. Отметим, что при наличии в моделях свободного газа и больших градиентах по давлению, эффективность метода АСРЯ падает, так как при этом доля ячеек множества Ое, возрастает. В худшем случае Ое = О и предобуславливатель ACPR становится эквивалентен предобуславливателю CPR.

Отличительной особенностью метода ACPR является его способность адаптироваться под условия в течении всего процесса численного решения.

а

СПИСОК ЛИТЕРАТУРЫ

1. Aziz, K. Pertoleum reservoir simulation / K. Aziz, A. Settari. Blitzprint Ltd., 2002. 480 p.

2. Wallis, J. R. Constraint residual acceleration of conjugate residual Method / J. R. Wallis, R. P. Kendall, T. E. Little // SPE 13536, 1985.

3. Stuben, K. Algebraic multigrid (AMG): An introduction with applications / K. Stuben // GDM Report 53, 03.1999.

4. Saad, Y. Iterative methods for sparse linear systems / Y. Saad, 2002.

5. Behie, A. Comparison of nested factorization, constrained pressure residual, and incomplete factorization preconditionings / A. Behie // SPE 13531-MS, 1985.

6. Coats, K. H. A note on IMPES and some IMPES based simulation models / K. H. Coats // SPE 65092, 2000.

ОБ АВТОРЕ

Борщук Олег Сергеевич, асп. каф. математики. Дипл. магистр прикл. математики и информатики (УГАТУ, 2005). Иссл. в обл. числ. моделир-я течения вязких, сжимаемых многофазн. жидкостей в пористой среде.