Численное решение задачи двухфазной фильтрации с неоднородными коэффициентами методом конечных элементов Текст научной статьи по специальности «Математика»

Математические заметки СВФУ Апрель—июнь, 2017. Том 24, № 2

УДК 519.632

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ДВУХФАЗНОЙ ФИЛЬТРАЦИИ С НЕОДНОРОДНЫМИ КОЭФФИЦИЕНТАМИ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ М. В. Васильева, Г. А. Прокопьев

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

Б01 10.25587/8УРи.2017.2.9245

Ключевые слова: пористая среда, двухфазная фильтрация, метод конечных элементов, метод Галеркина, численное моделирование.

Введение

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

Работа выполнена за счет Российского научного фонда, грант № 17—71—20055. © 2017 Васильева М. В., Прокопьев Г. А.

Основные уравнения фильтрации многофазной жидкости описываются законом сохранения массы и законом Дарси [4,8-11]. При рассмотрении модели многофазной фильтрации необходимо ввести такие понятия, как насыщенность Si и относительная проницаемость . Насыщенность ¿-й фазой флюида определяется как отношение объема, занимаемого данной фазой, к общему объему пустот в пористой среде. Будем рассматривать полностью насыщенную пористую среду:

^ 1.

5>

где Si — насыщенность ¿-й фазы флюида.

Закон сохранения массы для многофазной фильтрации описывается следующим дифференциальным оператором:

д(ФРгЗг) , ,

-—--Ь V • {ргЩ) = РгЧг,

где V — оператор дивергенции, ф — пористость пористой среды, ^ — плотность ¿-й фазы флюида, и = (и1,и2,и33) — скорость фильтрации ¿-й фазы флюида и ф — внешние источники или стоки. Отметим, что ^ отрицателен для стока и положителен для источников.

В дополнение к закону сохранения массы используем закон сохранения импульса в форме закона Дарси. Этот закон устанавливает линейное отношение между скоростью фильтрации и градиентом давления:

к

щ = ——кУрг,

у

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

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

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

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

Для численного решения возникающих задач для систем дифференциальных уравнений в частных производных используется смешанный метод конечных элементов, позволяющий напрямую выполнять законы сохранения [12]. Математическая модель записывается в смешанной обобщенной постановке в терминах скорость — давление и насыщенность. Для аппроксимации полей скоростей и давления могут быть использованы стандартные типы элементов, аналогичные рассматриваемым в [13-15] и др.

Уравнение насыщенности описывается уравнением конвекции-диффузии в случае наличия капиллярных сил и уравнением переноса при их отсутствии. Для аппроксимации уравнений конвекции-диффузии использование стандартного метода конечных элементов (метод Галеркина) ведет к возникновению ос-цилляций в решении задачи [16-18]. В отличие от оператора диффузионного переноса, аппроксимация которого ведет к симметричной и положительно определенной матрице, аппроксимация оператора конвекции с использованием стандартного метода Галеркина приводит к несимметричной знаконеопределенной матрице, что и приводит к возникновению осцилляций в решении задачи в случае задач с доминирующей конвекцией при больших числах Пекле. Стандартным способом борьбы с осцилляциями является измельчение расчетной сетки до тех пор, пока конвективный член не перестанет доминировать на сеточном уровне [5,19]. Отметим, что подобное измельчение сетки ведет к увеличению размерности разностной задачи, особенно для случая многомерной постановки, что также приводит к необходимости расчета задачи для скоростей и давления на соответствующих мелких сетках. Известно, что для борьбы с осцилляци-ями в методе конечных разностей используются потоковые схемы. Их можно построить посредством добавления искусственной диффузии для центрально-разностной схемы. Для конечно-элементной аппроксимации потоковые схемы можно построить аналогичным способом.

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

Рассмотрим течение двухфазной жидкости в пористой среде. Будем предполагать что пористая среда полностью насыщена двумя фазами:

1. Математическая модель двухфазной фильтрации

(1.1)

где Бы и Б о — насыщенность водой и нефтью соответственно.

'о

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

Рс = Ро - Р-ш, (1.2)

где ро и рш — давление нефтяной и водной фаз.

Учитывая отсутствие переноса массы между фазами (несмешиваемые жидкости), запишем закон сохранения массы для нефти и воды

д(фpoSo)

dt

д (фрш Sw)

+ V • (p0u0) = p0q0, + V • (pwUw) = Pw4w,

(1.3)

dt

где индексы o и w относятся к нефти и воде.

Для скорости течения имеем закон Дарси для каждой из фаз:

k k Uo = ——kVpo, uw = ——kVpw, (1.4)

po Pw

где kro = kro(So), krw = krw(Sw) и po, pw — эффективная проницаемость и вязкость для нефти и воды соответственно.

Для вывода математической модели введем понятие фазовой подвижности Ao и Aw

V k

j _ 'Ч д _ ^ri

А ' 1 рг'

где A = Ao + Aw — общая подвижность и fi — доля фазы в двухфазном потоке. Определим суммарную скорость двухфазного потока

U = Uo + Uw.

Сложив уравнения (1.3), с учетом соотношения (1.1) получим

-V- u = qo + qw, (1.5)

где для простоты предполагаем, что плотность воды и нефти постоянны, т. е. жидкость несжимаема, а также будем рассматривать несжимаемую пористую среду (ф = const).

Определим глобальное давление по формуле

p = Po + J fw(p'c(0) dt

Pc(So,Sw)

тогда

AVp = Aw Vpw + A oVpo.

Запишем фазовые скорости с использованием глобального давления:

uw = -k\w i Vp- yVPcj , (1.6)

u0 =-kX0 (vP - ^VPcy (1.7)

В качестве расчетной насыщенности будем использовать насыщенность водной фазой S = Sw, для которой имеем уравнение dS

Ф+ V(fwu) + V(kXoPyS) = qw, (1.8)

dt

гДеРс = it-

Для коэффициентов предположим следующие зависимости: qi = 0, kri = kri(Si), pi = const, = const, ф = const, i = o,w, и определим неоднородный тензор проницаемости пористой среды

к = k(x).

Таким образом, в некоторой области О математическая модель, описывающая процесс двухфазной фильтрации, задается системой уравнений для насыщенности, давления и скорости dS

ф— +u-(fi(S)VS)+V(kXop'c(S)VS) = 0, хеп, (1.9)

V^ u = 0, x e О, (1.10)

u = -k\(S)Vp, x e О, (1.11)

где US) = d-t-

Без учета капиллярных сил система уравнений принимает следующий упрощенный вид:

dS

+«.(/;(s-)vso = о, х&п, (1.12)

V- u = 0, x e о, (1.13)

u = -kX(S)Vp, x e О. (1.14)

Система уравнений дополняется начальным условием для насыщенности:

S(x, 0) = So, x e О, (1.15)

и граничными условиями для давления

dp

р = р 1, х е 1\, р = Р2, х е г2, —— = о, х е 90/Г1/Г2, (1-16)

дп

и насыщенности

= а:егь ^ = жег2. (1.17)

dx dx

Далее будем предполагать, что О — прямоугольная область, а 1\ и Г2 — часть

левой и правой границы соответственно (рис. 1).

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

Рис. 1. Расчетная область О

Рис. 2. Проницаемость пористой среды

Рис. 3. Графики зависимости капиллярных сил рс^) и их производных рС(£)

Рис. 4. Графики зависимости доли фазы воды fw ($) и ее производной frw ($)

2. Аппроксимация по времени

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

¿0 = 0, ¿„+1 = + т, п = 0,1, 2,....

В своем рассмотрении ограничимся двухслойными разностными аппроксимациями, которые связывают решение на п +1 временном слое с решением на п слое.

При решении нелинейных нестационарных задач математической физики естественно использовать линеаризованные схемы, в которых нелинейные коэффициенты берутся с предыдущего временного слоя. В этом случае адаптация к резким изменениям решения достигается выбором временных шагов. Неявная линеаризованная схема имеет следующий вид: С„+1 _ С„

ф--— +ип+1\7(,С(8п)\78п+1) + \7(к\оР'с(8п)\78п+1) = 0, (2.1)

т

V- и„+1 = 0, (2.2)

к-1ип+1 + Х(Бп)Урп+1 = 0. (2.3)

Полученные уравнения можно решать последовательно. Сначала рассчитываются давление и скорости (и, р), а затем — значения насыщенности Б.

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

• Улучшение устойчивости метода за счет неявного учета насыщенностей, но без совместного решения уравнений;

• Все функции, зависящие от насыщенности, к^(Б),рс(Б), вычисляются с использованием значений с предыдущего временного слоя.

3. Аппроксимация по пространству методом конечных элементов

Для аппроксимации по пространственным переменным применим метод конечных элементов и получим следующую вариационную формулировку: найти Б е ((О) и (и,р) е V(О) = [Я(О)]^ х £(О) такие, что

/с п+1 _ б п г

ф---гг1х+ Иип+\¡^{Бп)У Бп+1)г г1х

(кХ0р'с(Бп)ЧБп+1, Уг) йх = 0, (3.1)

о

У^ ип+1дйх = 0, (3.2)

J ((кА(Бп))-1 ип+1, V) йх _ ^ рп+1У • V йх + !рп+1 • ж йв = 0, (3.3)

о

?п))-1 ип+1 V йх _ I Рп+1У • v йх + j Рп+1

до

где г е ((О) и (V, д) е у (О) = [Я (О)]^ х Ь( О), й = 2, 3.

Для численного решения мы должны перейти от непрерывной вариационной задачи к дискретной для некоторого треугольного разбиения ¿РЬ области О, где ¡г — размерность расчетной сетки.

Введем конечномерные пространства (ъ С д < (ъ С ( и Vh С V, уъ С у и определим в них следующую дискретную задачу: найти Бъ е (ь и (иь ,рь) е такие, что

/сп+1 _ сп Г

_1 (кХор'^Б^УБ^1, Уг) йх = 0, (3.4)

У^ ип+1 д йх = 0, (3.5)

-9.916е+02 0 1172 2344 3.696еч-03

Рис. 5. Скорость фильтрации и = (их ,иу)

I((йА^))-1«)^1, V) ¿ж V ¿ж + 1 р1+1 • ™ ¿з = 0, (3.6)

О О до

где г е и («,д) е 14.

В качестве базисных функций для скорости и давления будем использовать элементы Равиа — Тома, а для насыщенности будем использовать линейные полиномы с дополнительной численной стабилизацией. Поскольку метод Галеркина на дискретном уровне дает аналог центрального разностного отношения конвективного члена уравнения, это ведет к возникновению осцилляций за счет немонотонности дискретного оператора конвекции. Известно, что для борьбы с осцилляциями в методе конечных разностей используются потоковые схемы. Противопотоковые схемы можно построить посредством добавления искусственной диффузии для центрально-разностной схемы.

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

Приведем результаты численного решения модельной задачи двухфазной фильтрации в прямоугольной области О (см. рис. 1). Проницаемость пористой среды зададим неоднородной из [20] (50 слой, рис. 2) и положим

Мо = 0.0003, м-ш = 0.003, ф = 1.

Рассмотрим три случая задания коэффициентов:

Случай 1. Линейные коэффициенты относительной проницаемости без учета капиллярных сил

ко = (1 - 5), к- =

1 = 1

2.000е-01 0.35 0,5 0.65 8.000е-01 И I I I I I ........... I

1 = 2.5

1 = 5

Рис. 6. Насыщенность 5 при линейном случае относительных проницаемостей (случай 1) на различных временных слоях

Случай 2. Нелинейные коэффициенты относительной проницаемости без

учета капиллярных сил

ко = (1 _ Б)2, кш = Б2

Случай 3. Нелинейные коэффициенты относительной проницаемости с учетом капиллярных сил

ко = (1 - Б)'2, кш = Б2, рс = РсьБ^ъ,

где в = 2.0 и рсь = 0.0001.

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

Б = 1, р = 1, х е гь Б = 0, р = 0 х е Г2.

Рис. 7. Насыщенность 5 при нелинейном случае относительных проницаемостей (случай 2) на различных временных слоях

Расчет проводился при Ттах = 0.20 с шагом по времени т = 0.004 с использованием 50 временных слоев. Численное моделирование проводились на неструктурированной треугольной равномерной расчетной сетке, содержащей 10181 вершин, 29937 ребер и 19958 ячеек.

Для проведения численных расчетов использовалась свободная библиотека ЕЕшС8 [21]. Для построения расчетной сетки использовалась свободно распространяемая программа СшвЬ [22]. Визуализация данных проводилась с использованием программы Paraview [23]. Для численного нахождения полей давления и скорости использовался прямой метод, который также был использован и при решении уравнения для насыщенности.

Графики зависимостей доли фазы воды и капиллярных сил и их производ-

Рис. 8. Численное исследование относительной погрешности для насыщенности в зависимости от изменения шага по времени т = 0.004, 0.008 и 0.01

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

Проведем численное сравнение зависимости решения задачи при различных значениях шага по времени т = 0.004, 0.008 и 0.01 (рис. 8). Мы вычисляли относительную погрешность

/(ие — и)2 ¿х п

е =

/и2 ¿х п

где ие — эталонное решение при т = 0.001. Сравнение проводилось для случая 3. Отметим, что погрешность решения при т = 0.004 менее десяти процентов, что оптимально с вычислительной точки зрения.

Распределения насыщенностей для нелинейных коэффициентов с учетом капиллярного давления представлены на рис. 9. При сравнении со случаем без учета капиллярных сил (рис. 10) получили более диффузирующее решение.

Ь = 1

Рис. 9. Насыщенность 5 при нелинейном случае относительных проницаемостей и наличии капиллярных сил (случай 3) на различных временных слоях

Рис. 10. Численное сравнение учета наличия капиллярных сил. Насыщенности для случая 1 (красное) и случая 2 (голубое) на последний момент времени (изолиния 5 = 0.25)

5. Заключение

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

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

ЛИТЕРАТУРА

1. Aziz K., Settari A. Petroleum reservoir simulation. Applied Sci. Publ., 1979.

2. Bear J. Dynamics of fluids in porous media. Dover Publ., 1988.

3. Chen Z., Huan G., Ma Y. Computational methods for multiphase flows in porous media. Society for Industrial Mathematics, 2006.

4. Афанасьева Н. М., Васильева М. В., Захаров П. Е. Параллельное численное моделирование процесса заводнения нефтяного месторождения // Мат. заметки СВФУ. 2011. Т. 18, вып. 1. С. 159-172.

5. Васильева М. В., Васильев В. И., Тимофеева Т. С. Численное решение методом конечных элементов задач диффузионного и конвективного переноса в сильно гетерогенных пористых средах // Уч. зап. Казан. ун-та. Сер. Физико-математические науки. 2016. Т. 158, № 2. С. 243-261.

6. Васильев В. И., Васильева М. В., Никифоров Д. Я. Решение задач однофазной фильтрации методом конечных элементов на вычислительном кластере // Вестн. Северо-Восточн. федерал. ун-та им. М. К. Аммосова. 2016. № 6. С. 8-17.

7. Akkutlu I. Y., Efendiev Y., Vasilyeva M., Wang Y. Multiscale model reduction for transport and flow problems in perforated domains. Multiscale model reduction for shale gas transport in poroelastic fractured media // J. Comput. Physics. 2018. V 353. P. 356-376.

8. Васильев В. И., Попов В. В., Тимофеева Т. С. Вычислительные методы в разработке месторождений нефти и газа. Новосибирск: Изд-во СО РАН, 2000.

9. Вабищевич П. Н. Явно-неявные вычислительные алгоритмы для задач многофазной фильтрации // Мат. моделирование. 2010. Т. 22, № 4. С. 118-128.

10. Vabishchevich P. N., Vasil'eva M. Iterative solution of the pressure problem for the multiphase filtration // Math. Modell. Anal. 2012. V. 17, N 4. P. 532-549.

11. Васильева М. В. Численное моделирование фильтрации на многопроцессорных системах // Мат. заметки ЯГУ. 2010. Т. 17, вып. 2. С. 105-112.

12. Афанасьева Н. М., Васильева М. В., Захаров П. Е. Параллельное численное моделирование процесса заводнения нефтяного месторождения // Мат. заметки ЯГУ. 2011. Т. 18,

вып. 1. С. 159-172.

13. Chung E. T., Leung W. T., Vasilyeva M. Mixed GMsFEM for second order elliptic problem in perforated domains // J. Comput. Appl. Math. 2016. V. 304. P. 84-99.

14. Brezzi F., Fortin M. Mixed and hybrid finite element methods. Berlin: Springer-Verl., 1991.

15. Raviart P. A., Thomas J. M. A mixed finite element method for 2-nd order elliptic problems // Mathematical aspects of finite element methods. Berlin; Heidelberg: Springer-Verl., 1977. P. 292-315.

16. Carstensen C. A posteriori error estimate for the mixed finite element method // Math. Comput. Amer. Math. Soc. 1997. V. 66, N 218. С. 465-476.

17. Самарский А. А., Вабищевич П. Н. Численные методы решения задач конвекции-диффузии. М.: Едиториал УРСС, 2004.

18. Donea J. Huerta A. Finite element methods for flow problems. Chichester: John Wiley & Sons Ltd, 2003.

19. Brooks A. N. A Petrov-Galerkin finite element formulation for convection dominated flows: PhD Thes. California Institute of Technology, 1981.

20. Алексеев В. Н., Васильева М. В., Степанов С. П. Итерационные методы решения для задачи течения и переноса в перфорированных областях // Вестн. Северо-Восточн. федерал. ун-та им. М. К. Аммосова. 2016. № 5. С. 67-79.

21. Christie M. A., M. J. Blunt. Tenth SPE comparative solution project: A comparison of upscaling techniques / SPE Reservoir Simulation Symposium. Society of Petroleum Engineers, 2001.

22. Logg A., Mardal K.-A., Wells G. N. Automated solution of differential equations by the finite element method. The FEniCS Book. 2011.

23. Software GMSH. (http://geuz.org/gmsh/).

24. Software package PARAVIEW. (http://www.paraview.org/).

Статья поступила 20 марта 2017 г.

Васильева Мария Васильевна, Прокопьев Григорий Анатольевич Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики, ул. Кулаковского 42, Якутск 677891 vasilyevadotmdotv@gmail.com, reilroot@gmail•com

Математические заметки СВФУ Апрель—июнь, 2017. Том 24, № 2

UDC 519.6

NUMERICAL SOLUTION TO THE PROBLEM OF TWO-PHASE FILTRATION WITH HETEROGENEOUS COEFFICIENTS

BY THE FINITE ELEMENT METHOD M. V. Vasilyeva and G. A. Prokopiev

Abstract. We consider the process of filtration of a two-phase fluid in a porous, heterogeneous medium. This process is described by a coupled system of equations for saturation, filtration rate, and pore pressure. We consider mathematical models with and without capillary forces, in the presence of which, for saturation, we have a nonstationary convection-diffusion equation. Since this process is characterized by a significant predominance of the convective term in the equation for saturation, counter-current approximations are used by adding non-uniform artificial diffusion. Speed and pressure are approximated using a mixed finite element method. The results of numerical calculations for a two-dimensional case with strongly heterogeneous permeability coefficients of a porous medium are presented. Several cases of relative fluid permeability associated with linear and nonlinear coefficients and the presence of capillary forces are considered.

DOI 10.25587/SVFU.2017.2.9245

Keywords: porous medium, two fasefiltration, finite elements method, Galerkin method, numerical simulation.

REFERENCES

1. Aziz K. and Settari A., Petroleum Reservoir Simulation, Appl. Sci. Publ. (1979).

2. Bear J., Dynamics of Fluids in Porous Media, Dover Publ. (1988).

3. Chen Z., Huan G., and Ma Y., Computational Methods for Multiphase Flows in Porous Media, Soc. Ind. Math. (2006).

4. Vasil'ev V. I., Popov V. V., and Timofeeva T. S., Vychislitel'nye Metody v Razrabotke Mestorozhdenij Nefti i Gaza [in Russian], Izdat. SO RAN, Novosibirsk (2000).

5. Vabishhevich P. N., Explicit-implicit computational algorithms for multiphase filtration problems," Math. Models Comput. Simul., 22, No. 4, 118-128 (2010).

6. Vabishhevich P. N. and Vasilyeva M. V., "Iterative solution to the pressure problem for the multiphase filtration," Math. Modelling Anal., 17, No. 4, 532-549 (2012).

7. Vasilyeva M. V., "Numerical modelling of filtration at multiprocessor systems," Mat. Zametki YaGU, 17, No. 2, 105-112 (2010).

8. Afanas'eva N. M., Vasil'eva M. V., and Zaharov P. E., "Parallel'noe chislennoe modelirovanie processa zavodnenija neftjanogo mestorozhdenija," Mat. Zametki YaGU, 18, No. 1, 159-172 (2011).

The work was supported by the grant from the Russian Scientific Foundation (code 17-7120055).

© 2017 M. V. Vasilyeva and G. A. Prokopiev

9. Chung E. T., Leung W. T., and Vasilyeva M., "Mixed GMsFEM for second order elliptic problem in perforated domains," J. Comput. Appl. Math., 304, 84-99 (2016).

10. Akkutlu I. Y., Efendiev Y., Vasilyeva M., and Wang Y., "Multiscale model reduction for transport and flow problems in perforated domains," J. Comput. Phys., 353, 356-376 (2018).

11. Vasilyeva M. V., Vasilyev V. I., and Timofeeva T. S., "Numerical solution by the finite elements method to the problems of transfer by diffusion and convection in strongly heterogenuous media," Uchen. Zap. Kazan. Univ., Ser. Fiz.-Mat. Nauki, 158, No. 2, 243-261 (2016).

12. Alekseev V. N., Vasilyeva M. V., and Stepanov S. P., "Iterative methods for solving the problems of flow and transfer in perforated domains," Vestn. Severo-Vostoch. Federal. Univ., No. 5, 67-79 (2016).

13. Vasil'ev V. I., Vasilyeva M. V., and Nikiforov D. Ya., "Solution to the problems of one-phase filtration by the finite elements method on the computation cluster," Vestn. Severo-Vostoch. Federal. Univ., No. 6 (2016).

14. Brezzi F. and Fortin M., Mixed and Hybrid Finite Element Methods, Springer, Berlin (1991).

15. Raviart P. A. and Thomas J. M., "A mixed finite element method for second order elliptic problems," in: Mathematical Aspects of Finite Element Methods, Springer, Berlin, Heidelberg, pp. 292-315 (1977).

16. Carstensen C. "A posteriori error estimate for the mixed finite element method," Math. Comput. Amer. Math. Soc., 66, No. 218, 465-476 (1997).

17. Samarskij A. A. and Vabishhevich P. N. Methods for Convection-Diffusion Problems [in Russian], Editorial URSS, Moscow (2004).

18. Donea J. and Huerta A., Finite Element Methods for Flow Problems, John Wiley & Sons Ltd, Chichester (2003).

19. Brooks A. N., A Petrov-Galerkin Finite Element Formulation for Convection Dominated Flows, Calif. Inst. Technology (1981).

20. Christie M. A. and Blunt M. J., Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques, SPE Reservoir Simulation Symp., Soc. Petroleum Engineers (2001).

21. Logg A., Mardal K.-A., and Wells G. N. (eds.), Automated Solution of Differential Equations by the Finite Element Method. The FEniCS Book, Springer, Heidelberg (2011).

22. Software GMSH. (http://geuz.org/gmsh/)

23. Software package PARAVIEW. (http://www.paraview.org/)

Submitted March 20, 2017

Maria V. Vasilyeva and Grigorii A. Prokopiev M. K. Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 42, Kulakovsky St., Yakutsk 677000, Russia vasilyeva_mv@mail .ru, reilroot@gmail .com