CC BY
8
2021
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Управление, вычислительная техника и информатика
№ 57
УДК 004.94 + 517.95 DOI: 10.17223/19988605/57/4
М.Г. Персова, Ю.Г. Соловейчик, И.И. Патрушев, А.С. Овчинникова
ПРИМЕНЕНИЕ ПРОЦЕДУРЫ ГРУППИРОВАНИЯ КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ ПОВЫШЕНИЯ ЭФФЕКТИВНОСТИ МОДЕЛИРОВАНИЯ НЕСТАЦИОНАРНОГО МНОГОФАЗНОГО ПОТОКА В ВЫСОКОНЕОДНОРОДНЫХ ТРЕХМЕРНЫХ
ПОРИСТЫХ СРЕДАХ
Работа выполнена при финансовой поддержке Министерства науки и высшего образования РФ в рамках проекта FSUN-2020-0012 (разработка общей вычислительной схемы моделирования трехмерного многофазного потока) и РФФИ в рамках научного проекта № 20-31-90049 (построение и исследование процедуры группирования).
Рассматривается проблема численного моделирования трехмерных многофазных нестационарных потоков в пористых средах применительно к задачам добычи углеводородов. Предлагается подход с неявным расчетом давления и явным переносом фаз и пересчетом насыщенностей на конечных элементах. Для повышения эффективности конечные элементы разбиваются на группы, в которых производится пересчет их состояния с допустимым и наиболее близким к требуемому шагом по времени.
Ключевые слова: многофазная фильтрация в пористых средах; моделирование нефтегазовых месторождений; метод конечных элементов.
Современные программные комплексы численного 3,0-моделирования процессов многофазной фильтрации в пористой среде широко используются при разработке нефтегазовых месторождений. Сложное строение коллекторов требует создания эффективных вычислительных схем, позволяющих проводить расчеты для моделей реальных месторождений с большим числом слоев, сильной неоднородностью среды и с большим количеством работающих скважин и зон перфораций. Успешность, а иногда и сама возможность решения таких важных для практики задач, как построение гидродинамических моделей месторождений путем решения обратных задач [1, 2] и синтез оптимального управления разработкой [3, 4], во многом зависит от эффективности вычислительных схем, применяемых при решении прямых задач.
Часто для решения задач фильтрации используют конечно-разностные численные схемы. В частности, подобные подходы используются в широко распространенных коммерческих симуляторах нефтедобычи, таких как Eclipse, Tempest и др. Однако многие авторы отмечают, что методы конечных разностей и конечных объемов обладают недостаточной геометрической гибкостью и недостаточной точностью при моделировании многофазных течений в высоконеоднородных средах [5, 6]. Публикуется довольно много работ, в которых для решения задач в сложных высоконеоднородных средах предлагается использовать различные модификации метода конечных элементов (МКЭ), например совместное использованием МКЭ и метода конечных объемов (FEFVM) [7, 8], mixed FEM и его специальные модификации (mixed hybrid FEM) [9-11]. Однако при использовании таких методов вычислительные затраты при решении ориентированных на практику задач нефтедобычи со сложной (многослойной латерально неоднородной) структурой среды и большим числом действующих скважин могут быть очень велики. Проблема снижения вычислительных затрат особенно актуальна при решении обратных задач, когда в процессе восстановления модели сложного месторождения требуется решать большое число соответствующих прямых задач.
В данной работе будет использован вариант МКЭ с непрерывными базисными функциями Гильбертова пространства функций с производными, суммируемыми с квадратом (Continuous Ga-
ЬгЫп [2, 12]). Необходимая точность выполнения законов сохранения масс достигается за счет специальной процедуры балансировки потоков [13]. Данная вычислительная схема была верифицирована на тесте SPE-10 [14] и применялась при моделировании реальных месторождений высоковязкой нефти Республики Татарстан [2]. Повышение вычислительной эффективности для задач такого класса достигается в том числе и за счет использования специальных некомфорных сеток [15], которые позволяют существенно сократить число степеней свободы в конечноэлементных аппроксимациях без увеличения погрешности численного решения. Свою эффективность предложенный подход показал при решении обратных задач [2] и при синтезе оптимального управления разработкой месторождения [4].
Принцип построения предлагаемой вычислительной схемы аналогичен принципу известных схем IMPES [16] и 1МРЕС [10, 11] с конечноэлементным (неявным) расчетом давления и явным переносом фаз между ячейками конечноэлементной сетки на каждом временном шаге. Вычислительные затраты напрямую зависят от величины временного шага, дробление которого может потребоваться при большой неоднородности среды и сильно изменяющейся в пространстве скорости потока.
В данной работе будут представлены схема группирования конечных элементов по временному шагу и алгоритм моделирования перетоков фаз с учетом разбиения ячеек сетки на группы, позволяющие минимизировать количество вычислений поля давления, не понижая качества аппроксимации многофазного потока. В качестве теста использовалась задача сравнительного проекта SPE-10. Показана сходимость данного метода, проведены сравнение с данными участников проекта SPE-10, анализ затрат машинного времени и точности получаемых решений при разных временных шагах.
1. Математическая модель
В задачах моделирования месторождений расчет движения многофазной смеси выполняется в неоднородной пористой среде, которая характеризуется зависящими от пространственных координат тензором структурной проницаемости К и пористостью Ф. Скорость ут движения в пористой среде каждой фазы фильтрующейся смеси подчиняется закону Дарси:
т { т Л
где т - номер фазы, V™ - скорость потока фазы, Р - давление, Р'" - капиллярное давление, g -ускорение свободного падения, кт - коэффициент относительной фазовой проницаемости, л™ - динамическая вязкость, рт - плотность фазы.
Распределение давления в расчетной области О описывается краевой задачей [2]:
^ ИР кт { / \ / \ТЛ^ ИР
= !/т,°, (1)
I Кт К ( grad (Р + Рст) + (0,0, рmg)
и=1 л V ' 4 '
т=1
I К™ К ( grad (Р + Рст) + (0,0, р^)
п=1 л V
ИР „
■п=^Г'Т. (3)
т=1
т=1 л
Р|Г1 = РГ, (2)
ИР
I
т=1 Л
Здесь ИР - количество фаз. Функции ут'° соответствуют объемным источникам (стокам) фаз в расчетной области О, возникающим, например, вследствие выделения / поглощения газа жидкими фазами, химических процессов или сжатия / расширения фаз при изменении давления. Г1 - объединение тех границ расчетной области О, где задано давление РГ, а Г2 - объединение границ О, где
Г ИР Г Г 2
задан поток смеси / =1 /т, . Функция / не равна нулю на тех границах из Г , которые соот-
т=1
ветствуют активным (в интервале времени At) зонам перфорации. Остальные границы из Г2 являются непроницаемыми (на них /Г = 0).
Краевая задача (1)-(3) решается методом конечных элементов на шестигранных неконформных сетках [15]. По полученным значениям давления вычисляются объемы смеси, перетекающие через грани Гг конечных элементов Ое за единицу времени:
г,.
'е
№ кт ( / \ / \Т
I ^т к( (Р+рт)+(о,о, ?тё)
п=1Л V
?Г о ¿г,
т=1 Л
где йр ,п ~ внешняя (по отношению к ) нормаль к Гг. Для внутренних граней Гг объем смеси, перетекающий за единицу времени из конечного элемента (ячейки) Ое в смежную с ним по грани Гг ячейку 0.к (Гг = О, П 0.к) определяется как взвешенное среднее объемов, перетекающих через эту границу на конечных элементах Ое и Ок :
0Г = Хк ОГ о + Хе ОГ о . (4)
_ т _
Коэффициенты А, в (4) ('/,к для О,, и /.„ для С1е) определяются как /. = ^^~7• гДе ^ = "г ■
ш=1 Л
т.е. /С определяется по значению тензора К на соответствующем конечном элементе.
Если для решения краевой задачи (1)-(3) применяется МКЭ с базисными функциями из Н1 (так называемый CG [12]), то получается численное решение, не гарантирующее сохранения масс веществ в фильтрующейся смеси [7, 17] (закон сохранения в этом случае лишь аппроксимируется с той или иной точностью в зависимости от подробности сетки). Поэтому мы используем специальный метод балансировки потоков [13], который корректирует перетекающие объемы ОГ таким образом, чтобы законы сохранения масс отдельных фаз (и их компонент) были выполнены с необходимой точностью. По сбалансированным потокам смеси осуществляются перетоки фаз через Гг и вычисляется новый фазовый состав в конечных элементах. Для этого определяются объемы фаз От, перетекающих через грань Г; в единицу времени. В условиях, когда эффект гравитации или капиллярного давления является существенным, разные фазы могут перетекать через грань Г в противоположных направлениях. С учетом этого алгоритм выполнения перетоков между конечными элементами выглядит следующим образом.
Для всех конечных элементов Ое и принадлежащих им граней Г вычисляются численные потоки фаз От, о по формуле
т / -р \
От,,«. =-|—к|ёгас1(/^/Г) + ((),(),р'"я) (5)
Уравнение (5) определяет не только величину потока т-й фазы через грань Г , но и его направление по отношению к конечному элементу Ое. Фаза т вытекает из конечного элемента Ое, если значение Ог ~ & положительное, и втекает при отрицательном значении. Численный поток фазы От через грань Г, берется с подветренной стороны:
= 0,0 е- ОГ.Ое > 0. (6)
' I Ог,,о если от,о< 0.
Затем для граней Гг для каждой фазы вычисляются величины П^ , которые фактически определяют долю т-й фазы в потоке смеси через грань Г;:
ШР,
щ = о™ [/¿ог, |. (7)
Далее мы корректируем численные потоки фаз ОТ таким образом, чтобы сумма потоков всех фаз была строго равна сбалансированному потоку смеси ()г , перетекающему через грань Гг. Для этого вычислим разницу между численным и сбалансированным потоками смеси и распределим ее между потоками фаз пропорционально величинам ОТ , т.е. перевычислим 0Т по формуле
вг:=вг:+[вг,-хвг,}^. (8)
Таким образом, объем т-й фазы V™, который за время & перетекает через грань Гг, вычисляется с использованием по формуле
V™ = 0т, , (9)
и новые значения насыщенностей на каждом элементе Пе на конец шага по времени &
( \
ri m _
Sm mes(Ое)0 + AFQm + Z VTm - Z VTm /(mes(Ое)ф), (10)
^ ^ Tin ni ^ .-mit m ' / * '
iel™ ■ ieT,
где mes (Qe ) - объем ячейки fie ; AVm - дефицит / профицит объема m-й фазы в конечном элементе m, который может образоваться, например, в результате химических процессов или сжатия / расширения фаз при изменении давления; iout,m, i™'m - множества номеров граней элемента Qe, через которые m-я фаза вытекает из Qe и, соответственно, втекает в Qe.
На основе полученных значений S^ вычисляются новые значения фазовых проницаемостей
к" согласно заданным зависимостям к" от насыщенностей фаз. Затем осуществляется переход
к следующему шагу по времени, на котором процедура повторяется, начиная с расчета давления.
Теперь более подробно остановимся на выборе временного шага At. Кроме того, что он непосредственно влияет на точность аппроксимации по времени, его значение должно быть ограничено величиной объема фаз в ячейках, из которых эти фазы вытекают. Шаг At должен быть таким, чтобы в каждой ячейке суммарный вытекающий объем фазы m не превышал имеющийся объем подвижной фазы в ней. Это естественное условие определяет его предельное (максимальное) значение:
At <((S" - S"'res ) mes (m )Ф + A V" )/ Z Qm. , V^, Vm , (11)
\\e e ' e H j ey-out,m i
/ me
о m res i
где - остаточная насыщенность фазы m в ячейке ile.
Таким образом, временной шаг At зависит от размеров конечных элементов, пористости, величин Sm и потоков фаз Vm, перетекающих через грани конечных элементов, и для тех элементов,
me А i
размеры которых невелики и / или через которые перетекают большие потоки, шаг At может быть очень маленьким. И если такой шаг использовать для обработки всех ячеек, затраты машинного времени на осуществление перетоков фаз могут быть очень большими.
2. Группирование конечных элементов и алгоритм расчета перетоков фаз
В данной работе при выполнении процедуры перетока фаз между ячейками мы будем распределять ячейки конечноэлементной сетки по группам, в каждой из которых может быть использован свой временной шаг, удовлетворяющий критерию (11). Данная процедура позволит выбрать некоторый глобальный временной шаг Atmain, определяемый требованиями к качеству аппроксимации по времени и не зависящий от объемов фаз в отдельных ячейках.
Обозначим через множество номеров конечных элементов, определяющих группу ячеек О = , е е ^ |, для которых величина временного шага Д^ определяется соотношением Д^ = ДТшп/22—1 , т.е. шаг по времени для первой группы (О) равен Дгтш п, а шаг по времени каждой группы 02 для 2 > 1 в 22—1 раза меньше Дгтш п . Количество таких групп обозначим N° .
Для того, чтобы распределить ячейки Ое по группам, для каждой т-й фазы в ячейке Ое согласно условию (11) определяем временной шаг, являющийся допустимым для этой фазы:
(12)
^ а =(( sm -
Ф+Avm )/ s Qm.
iel 0U"m '
') mes (Qe )
Номер группы для конечного элемента Qe выбирается из условия
minjg : Aîg < min AtmA
(, g m
(13)
Перенос фаз между ячейками осуществляется следующим образом. Обработка ячеек начинается с группы с самым большим номером 2 = N° (т.е. с самым маленьким шагом по времени). Делается два шага по времени Д для всех ячеек этой группы. Затем делается один шаг по времени для ячеек
о
группы с номером 2 = N — 1, после чего опять делается два шага по времени для ячеек группы с ноС О "ГТ
мером N . Далее делается еще один шаг по времени для ячеек с номером группы 2 = N — 1. После этой второй обработки данной группы делается первый шаг обработки группы с номером 2 = N° — 2 и опять осуществляется возврат к обработке группы с номером 2 = N° .
Рис. 1. Алгоритм, реализующий расчет перетоков фаз Fig. 1. Algorithm of phase transport
Таким образом, после первой обработки каждой группы с номером g осуществляется возврат к обработке ячеек группы с номером N° (и далее до группы g), а после второй обработки группы с номером g осуществляется переход к группе g — 1 (с обнулением счетчиков СигО обработки групп с большими номерами). Алгоритм завершается после обработки группы с номером 1 (т.е. группы с самым большим шагом А/г ).Обработка ситуации, когда фазы перетекают между ячейками из
разных групп, осуществляется следующим образом. Прежде всего отметим, что при вычислении перетекающих объемов фаз VТ по формуле (9) всегда используется шаг Агг группы того элемента, из
которого вытекает поток 0т . При пересчете насыщенностей в ячейках это позволяет использовать вытекающие из них объемы фаз V™ без изменений. А втекающие объемы вычисляются в зависимо-
Г I
сти от того, какой группе принадлежит соседний элемент, из которого фазы перетекают в обрабатываемую ячейку. Таким образом, если при обработке некоторого конечного элемента из группы через одну из его граней втекала фаза из ячейки группы с меньшим или равным номером г , то перетекающий объем VТ делится на 2g—г, чтобы учесть разницу во временных шагах. Если же в элемент
Г I
из группы Gg втекает фаза из ячейки группы с большим номером г , то используется ранее накопленное значение перетекающего объема УштТ , которое было получено в ходе предыдущих обработок
ячейки группы с номером г .
Алгоритм, реализующий процедуру перетоков смеси для расчета нового состояния ячеек с учетом группирования, изображен на рис. 1.
3. Численные эксперименты
Продемонстрируем работоспособность и эффективность процедуры группирования на примере десятой тестовой задачи из сравнительного проекта SPE, рассмотренной в [14]. Модель представляет собой куб с размерами 1200 х 2200 х 170 фут3 (365,76 х 670,56 х 51,816 м3), в углах которого расположены четыре добывающие скважины (P1-P4), а в центре находится нагнетательная скважина. Добывающие скважины работают при заданном в них давлении 272 атмосферы, а нагнетательная скважина закачивает 5 000 баррелей воды в сутки. Исходная модель содержит 85 слоев и характеризуется высокой неоднородностью структурной проницаемости и пористости. Значения параметров задачи заданы на регулярной сетке, содержащей 1 122 000 ячеек (60 х 220 х 85 ячеек), и варьируют в диапазоне от 6,65 • 10-4 до 20 000 мД для коэффициентов тензора проницаемости и в диапазоне от 1,4 10-6 до 0 для пористости.
Для этой модели была выполнена предварительная обработка с укрупнением ячеек сетки по вертикали (upscaling). В результате расчет проводился на сетке из 244 052 ячеек. Время жизни «месторождения» - 2 000 суток.
Для оценки эффективности предлагаемого нами подхода с группированием ячеек была проведена серия расчетов с разными временными шагами Atmam, равными 100, 50, 20, 10, 5, и 1 сут. Также был проведен расчет без использования процедуры группирования. В нем расчет давления выполнялся с шагом А = 20 суток, а перетоки фаз - с одинаковым для всех ячеек временным шагом, удовлетворяющим условию (11). На рис. 2 представлены соответствующие результаты для скважин P1 и P3 (рис. 3): дебит нефти в баррелях (см. рис. 2, а) и обводненность (см. рис. 2, b) для шагов Atmai", равных 100, 50, 20 и 1 сут. Очевидно, с дроблением шага по времени наблюдается сходимость. Данные дебита нефти на скважине P3 при шаге Atmain на шагах 20, 10 и 5 сут отличаются от значений при Atmain=1 сутки в среднем на 2, 1 и 0,5% соответственно, а на скважине P1 - 1, 0,5 и 0,2%. Также от-
метим, что дебит нефти при расчете с шагом Агташ=1 сут очень хорошо совпадал с результатами, полученными без применения процедуры группирования. Аналогичная картина сходимости наблюдается и для обводненности добываемой смеси.
«33. 1333. 1500 fdv) 5-33. 10». ИМ г(itpy
a b
Рис. 2. Графики дебита нефти (а) и обводненности (b) для расчетов с различными Atmain : квадрат - 100 сут, незакрашенный круг - 50 суток, треугольник - 20 суток, закрашенный круг - 1 сутки Fig. 2. Oil rate (a) and water cut (b) curves for calculation with different Atmal" : square - 100 days, open circle - 50 days, triangle - 20 days, filled circle - 1 day
a) b)
Рис. 3. Распределение ячеек по группам (зеленый цвет - группа с большим номером, желтый цвет - с меньшим) Fig. 3. Distribution of cells into groups (green is the group with a higher number, yellow is the group with a lower number)
Временные затраты для расчетов с различным шагом Atma'n
Atmain, сут Количество групп Расчет давления Расчет потоков Расчет перетоков Общее время расчета
1 8 13 ч 56 мин 1 ч 55 мин 1 ч 14 мин 17 ч 37 мин
5 10 2 ч 40 мин 24 мин 31 мин 3 ч 42 мин
10 11 1 ч 25 мин 13 мин 29 мин 2 ч 12 мин
20 12 43 мин 7 мин 27 мин 1 ч 20 мин
50 14 17 мин 3 мин 26 мин 48 мин
100 15 8 мин 2 мин 26 мин 38 мин
Без группирования 43 мин 6 мин 66 ч 46 мин 84 ч 31 мин
В таблице приведены вычислительные затраты основных процедур описанной вычислительной схемы. На каждом временном шаге Atmain производится решение двух СЛАУ: первой - при решении краевой задачи (1)-(3) методом конечных элементов (размер СЛАУ определяется количеством узлов сетки), второй - в ходе выполнения процедуры балансировки для расчета потоков смеси (размер этой СЛАУ определяется количеством граней [13]). Для решения обоих СЛАУ применялся прямой решатель PARDISO из библиотеки Intel MKL.
Сопоставляя погрешность по времени с вычислительными затратами при различных значениях Atmain, можно считать оптимальным значение Atmam = 20 суток. Для этого расчета на рис. 3 приведено распределение ячеек по группам на конец расчета, где желтый цвет соответствует группе с наименьшим номером (шаг по времени AtG = 20 сут), а темно-зеленый цвет соответствует ячейкам из
группы с наибольшим номером (шаг по времени Air = 28 мин). Из рисунка видно, что ячейки, расположенные вдоль основных направлений течения (от нагнетательной скважины I1 к добывающим P1, P2, P3, P4) и в зонах повышенной проницаемости, попадают в группы с большими номерами, и основные вычислительные затраты при расчете перетоков приходятся на них. Если же не использовать процедуру группирования, то затраты машинного времени составляют порядка 85 час из-за условия Куранта-Фридрихса-Леви (CFL) и связанной с ним необходимости задания довольно мелкого шага по времени для выполнения процедуры перетоков фаз (At примерно 20-30 мин).
Заключение
Рассмотрена вычислительная схема моделирования многофазных потоков в пористых средах при решении задач нефтедобычи. Схема основана на неявном расчете давления с использованием МКЭ, балансировке численных потоков смеси и явном переносе фаз между конечными элементами. Предложена процедура группирования и основанный на ней алгоритм переноса фаз, позволяющие выполнять пересчет насыщенностей в ячейках с разными шагами по времени в зависимости от скорости течения и объема фаз в них. Тем самым устраняется необходимость использования слишком мелких шагов по времени из-за возможного нарушения условия CFL на отдельных ячейках, и резко снижаются вычислительные затраты без ущерба точности получаемого решения.
Эффективность предложенной процедуры группирования и алгоритма переноса фаз продемонстрирована на тесте SPE-10, в котором нефтедобыча моделируется в высоконеоднородной среде. Для этой задачи проведенные вычислительные эксперименты показали, что за счет использования предложенной процедуры группирования время счета сокращается примерно на полтора порядка без снижения точности получаемого решения.
ЛИТЕРАТУРА
1. Bukshtynov V., Volkov V., Durlofsky O., Aziz K. Comprehensive framework for gradient-based optimization in closed-loop
reservoir management // Comput. Geosci. 2015. V. 19, № 4. P. 877-897.
2. Persova M.G., Soloveichik Yu.G., Vagin D.V., Grif A.M., Kiselev D.S., Patrushev I.I., Nasybullin A.V., Ganiev B.G. The design
of high-viscosity oil reservoir model based on the inverse problem solution // J. Pet. Sci. Eng. 2021. V. 199. Art. 108245.
3. Ni H.M., Liu Y.J., Fan Y.C. Optimization of injection scheme to maximizing cumulative oil steam ratio based on improved
artificial bee colony algorithm // J. Pet. Sci. Eng. 2019. V. 173. P. 371-380.
4. Persova M.G., Soloveichik Yu.G., Vagin D.V., Grif A.M., Patrushev I.I., Ovchinnikova A.S. Oil production optimization based on
the finite-element simulation of the multi-phase flow in porous media and inverse problem solution // Irkutsk: GeoBaikal 2020 (EAGE). 2020. V. 2020, № 1. P. 1-6.
5. Doyle B., Riviere B., Sekachev M. A multinumerics scheme for incompressible two-phase flow // Comput. Methods Appl. Mech.
Eng. 2020. V. 370. Art. 113213.
6. Abd A.S., Abushaikha A. Velocity dependent up-winding scheme for node control volume finite element method for fluid flow in
porous media // Sci. Rep. 2020. V. 10. Art. 4427.
7. Schmid K.S., Geiger S., Sorbie K.S. Higher order FE-FV method on unstructured grids for transport and two-phase flow with
variable viscosity in heterogeneous porous media // J. Comput. Phys. 2013. V. 241. P. 416-444.
8. Abushaikha A.S., Blunt M.J., Gosselin O.R., Pain C.C., Jackson M.D. Interface control volume finite element method for
modelling multi-phase fluid flow in highly heterogeneous and fractured reservoirs // J. Comput. Phys. 2015. V. 298. P. 41-61.
9. Zhang N., Yan B., Sun Q., Wang Y. Improving multiscale mixed finite element method for flow simulation in highly heterogeneous
reservoir using adaptivity // J. Pet. Sci. Eng. 2017. V. 154. P. 382-388.
10. Moortgat J., Firoozabadi A. Higher-order compositional modeling of three-phase flow in 3D fractured porous media based on cross-flow equilibrium // J. Comput. Phys. 2013. V. 250. P 425-445.
11. Amooie M.A., Moortgat J. Higher-order black-oil and compositional modeling of multiphase compressible flow in porous media // Int. J. Multiph. Flow. 2018. V. 105. P. 45-59.
12. Scovazzi G., Wheeler M.F., Mikelic A., Lee S. Analytical and variational numerical methods for unstable miscible displacement flows in porous media // J. Comput. Phys. 2017. V. 335. P. 444-496.
13. Persova M.G., Soloveichik Yu.G., Grif A.M., Patrushev I.I. Flow Balancing in FEM Modelling of Multi-Phase Flow in Porous Media // Novosibirsk: 2018 XIV Int. Scientific-Technical Conference on Actual Problems of Electronics Instrument Engineering (APEIE). 2018. P. 205-211.
14. Christie M.A., Blunt M.J. Tenth SPE comparative solution project: A comparison of upscaling techniques // SPE Reservoir Evaluation and Engineering. 2001. V. 4, № 4. P. 308-316.
15. Persova M.G., Soloveichik Y.G., Vagin D.V., Kiselev D.S., Koshkina Yu.I. Finite element solution to 3-D airborne time-domain electromagnetic problems in complex geological media using non-conforming hexahedral meshes // J. Appl. Geophys. 2020. V. 172. Art. 103911.
16. Chen H., Kou J., Sun S., Zhang T. Fully mass-conservative IMPES schemes for incompressible two-phase flow in porous media // Comput. Methods Appl. Mech. Eng. 2019. V. 350. P. 641-663.
17. Zhang R.H., Zhang L.H., Luo J.X., Yang Z.D., Xu M.Y. Numerical simulation of water flooding in natural fractured reservoirs based on control volume finite element method // J. Pet. Sci. Eng. 2016. V. 146. P. 1211-1225.
Поступила в редакцию 21 мая 2021 г.
Persova M.G., Soloveichik Yu.G., Patrushev I.I., Ovchinnikova A.S. (2021) APPLICATION OF THE FINITE ELEMENT GROUPING PROCEDURE TO IMPROVE THE EFFICIENCY OF UNSTEADY MULTIPHASE FLOW SIMULATION IN HIGH-HETEROGENEOUS 3D POROUS MEDIA. Vestnik Tomskogo gosudarstvennogo universiteta. Upravlenie, vychislitelnaja tehnika i informatika [Tomsk State University Journal of Control and Computer Science]. 57. pp. 34-44
DOI: 10.17223/19988605/57/4
The problem of three-dimensional modeling of oil production in highly heterogeneous reservoirs is considered. A semi-implicit scheme similar to IMPES or IMPEC is used for modeling. The proposed approach is based on implicit calculation of pressure using finite element method (Continuous Galerkin with flow balancing), and explicit phase transport and recalculation of their saturations on finite elements. High computational efficiency is achieved due to the division of finite elements into groups, in which an explicit recalculation of their state (new saturation values after phase transport) is performed with its own (for each group) allowable time step closest to the required.
The cells fie of the finite element mesh are divided into groups as follows. Let Atmaln be the time step determined only by the requirements of the approximation quality and not depending on the volume of the phases flowing between the finite elements.
Let us denote by the set of numbers of finite elements defining a group Gg = jfie, e e Ig J of cells for which the value of the
time step AfG is determined by the ratio AtG = Atmam/2g, i.e. the time step of each group Gg is 2g_1 times less than the "global"
step Atmam . The number of groups will be denoted as NG .
We divide the finite elements into groups so that there is enough of each m-th phase in fie for its transport. To do this, for each finite element and all phases present in it, the time step is determined as:
A^p, =(fe "Sm,res)mes(Qe)0 + AVme)/ Z |Qrm|,
where Sp and Sm'res are the saturation and residual saturation of the m-th phase in the cell p , respectively; mes(p) is the volume of p ; ® is the porosity; AVpf is deficit/surplus of the m-th phase in the cell Pe, which appears, for example, due to gas dissolution/evolution, compression/expansion of phases caused by pressure changes, or as a result of chemical processes; I0ut'm is the set of numbers of element faces through which the m-th phase flows out of p , Qm is the flow of the m-th phase through the face r,-. The group number for the finite element Pe is chosen from the condition
rnnjg : At0 < min Atm p J .
\ s m )
Phase flows from one cell to another in accordance with step AiG increase. First the cells from the groups with larger numbers g
are processed. The transition to the group with a lower number g is carried out when such a number of steps have been passed in the groups with higher numbers that their sum has reached the values AiG . For this, after processing each group with the number g, the
transition is made either to the group with the number g - 1, if the accumulated time step has reached the value Ata i , or to the
group with the number g = NG .
The situations in which the phases flow between cells from different groups are processed as follows. First of all+ we note that when calculating the flowing phase volumes Vrm we use the group step AiG of the element the phase flows out from. This makes it
possible to use the phase volumes Vrm flow out of these cells without changes when recalculating the saturations in the cells. And the
inflowing volumes are calculated depending on the group number of the neighboring element the phases flow into the processed cell from.
Thus, when processing some finite element from the g-th group, if a phase inflows through one of its faces from a cell of the
group with a smaller or equal number r, then the flowing volume Vrm is divided by 2g~r to take into account the difference between
the time steps. On the other hand, if a phase flows into an element of the g-th group from the cell of the group with a larger number r,
then we use the previously accumulated value Vsumm of the flowing volume, which has been obtained as the result of the previous
processing of the cell from the group with the number r.
Computational experiments were carried out using the SPE-10 model. The obtained results indicated that for this problem the computational costs can be reduced by one and half orders due to the grouping of finite elements and the use of own allowable steps in groups.
Keywords: multiphase flow in porous media; modeling of oil and gas fields; finite element method.
PERSOVA Marina Gennadievna (Doctor of Technical Sciences, Professor, Novosibirsk State Technical University, Novosibirsk, Russian Federation). E-mail: persova@ami.nstu.ru
SOLOVEICHIK Yuri Grigorievich (Doctor of Technical Sciences, Professor, Novosibirsk State Technical University, Novosibirsk,
Russian Federation).
E-mail: soloveychik@ami.nstu.ru
PATRUSHEVIlya Igorevich (Post-graduate Student, Novosibirsk State Technical University, Novosibirsk, Russian Federation). E-mail: patrushev@ami.nstu.ru
OVCHINNIKOVA Anastasia Sergeevna (Post-graduate Student, Novosibirsk State Technical University, Novosibirsk, Russian Federation).
E-mail: ovchinnikova.2014@stud.nstu.ru
REFERENCES
1. Bukshtynov, V., Volkov, V., Durlofsky, O. & Aziz, K. (2015) Comprehensive framework for gradient-based optimization in
closed-loop reservoir management. Computer & Geosciences. 19(4). pp. 877-897. DOI: 10.1007/s10596-015-9496-5
2. Persova, M.G., Soloveichik, Yu.G., Vagin, D.V., Grif, A.M., Kiselev, D.S., Patrushev, I.I., Nasybullin, A.V. & Ganiev, B.G.
(2021) The design of high-viscosity oil reservoir model based on the inverse problem solution. Journal of Petroleum Science and Engineering. 199. Art. 108245. DOI: 10.1016/j.petrol.2020.108245
3. Ni, H.M., Liu, Y.J. & Fan, Y.C. (2019) Optimization of injection scheme to maximizing cumulative oil steam ratio based
on improved artificial bee colony algorithm. Journal of Petroleum Science and Engineering. 173. pp. 371-380. DOI: 10.1016/j.petrol.2018.10.032
4. Persova M.G., Soloveichik Yu.G., Vagin D.V., Grif A.M., Patrushev I.I. & Ovchinnikova A.S. (2020) Oil production optimization
based on the finite-element simulation of the multi-phase flow in porous media and inverse problem solution. Irkutsk: GeoBaikal 2020 (EAGE). 2020 (1). pp. 1-6. DOI: 10.3997/2214-4609.202052021
5. Doyle, B., Riviere, B. & Sekachev, M. (2020) A multinumerics scheme for incompressible two-phase flow. Computer Methods
in Applied Mechanics and Engineering. 370. Art.113213. DOI: 10.1016/j.cma.2020.113213
6. Abd, A.S. & Abushaikha, A. (2020) Velocity dependent up-winding scheme for node control volume finite element method
for fluid flow in porous media. Scientific Reports. 10. Art. 4427. DOI: 10.1038/s41598-020-61324-4
7. Schmid, K.S., Geiger, S. & Sorbie, K.S. (2013) Higher order FE-FV method on unstructured grids for transport and two-phase
flow with variable viscosity in heterogeneous porous media. Journal of Computational Physics. 241. pp. 416-444. DOI: 10.1016/j.jcp.2012.12.017
8. Abushaikha, A.S., Blunt, M.J., Gosselin, O.R., Pain, C.C. & Jackson, M.D. (2015) Interface control volume finite element method
for modelling multi-phase fluid flow in highly heterogeneous and fractured reservoirs. Journal of Computational Physics. 298. pp. 41-61. DOI: 10.1016/j.jcp.2015.05.024
9. Zhang, N., Yan, B., Sun, Q. & Wang, Y. (2017) Improving multiscale mixed finite element method for flow simulation in highly
heterogeneous reservoir using adaptivity. Journal of Petroleum Science and Engineering. 154. pp. 382-388. DOI: 10.1016/j.petrol.2017.04.012
10. Moortgat, J. & Firoozabadi, A. (2013) Higher-order compositional modeling of three-phase flow in 3D fractured porous media based on cross-flow equilibrium. Journal of Computational Physics. 250. pp. 425-445. DOI: 10.1016/j.jcp.2013.05.009
11. Amooie, M.A. & Moortgat, J. (2018) Higher-order black-oil and compositional modeling of multiphase compressible flow in porous media. International Journal of Multiphase Flow. 105. pp. 45-59. DOI: 10.1016/j.ijmultiphaseflow.2018.03.016
12. Scovazzi, G., Wheeler, M.F., Mikelic, A. & Lee, S. (2017) Analytical and variational numerical methods for unstable miscible displacement flows in porous media. Journal of Computational Physics. 335. pp. 444-496. DOI: 10.1016/j.jcp.2017.01.021
13. Persova, M.G., Soloveichik, Yu.G., Grif, A.M. & Patrushev, I.I. (2018) Flow Balancing in FEM Modelling of Multi-Phase Flow in Porous Media. Novosibirsk: 2018 XIV Int. Scientific-Technical Conference on Actual Problems of Electronics Instrument Engineering (APEIE). pp. 205-211. DOI: 10.1109/APEIE.2018.8545457
14. Christie, M.A. & Blunt, M.J. (2001) Tenth SPE comparative solution project: A comparison of upscaling techniques. SPE Reservoir Evaluation and Engineering. 4(4). pp. 308-316. DOI: 10.2118/66599-MS
15. Persova, M.G., Soloveichik, Y.G., Vagin, D.V., Kiselev, D.S. & Koshkina, Yu.I. (2020) Finite element solution to 3-D airborne time-domain electromagnetic problems in complex geological media using non-conforming hexahedral meshes. Journal of Applied Geophysics. 172. Art. 103911. DOI: 10.1016/j.jappgeo.2019.103911
16. Chen, H., Kou, J., Sun, S. & Zhang, T. (2019) Fully mass-conservative IMPES schemes for incompressible two-phase flow in porous media. Computer Methods in Applied Mechanics and Engineering. 350. pp. 641-663. DOI: 10.1016/j.cma.2019.03.023
17. Zhang, R.H., Zhang, L.H., Luo, J.X., Yang, Z.D. & Xu, M.Y. (2016) Numerical simulation of water flooding in natural fractured reservoirs based on control volume finite element method. Journal of Petroleum Science and Engineering. 146. pp. 1211-1225. DOI: 10.1016/j.petrol.2016.08.024
