CC BY
65
____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА
Том 155, кн. 3 Физико-математические науки
2013
УДК 532.546
РЕШЕНИЕ ЗАДАЧ ТРЕХФАЗНОЙ ФИЛЬТРАЦИИ ЖИДКОСТИ В ТРЕХМЕРНЫХ ПЛАСТАХ С ГИДРОДИНАМИЧЕСКИ НЕСОВЕРШЕННЫМИ СКВАЖИНАМИ НА ГЕТЕРОГЕННЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ
А.В. Цепаев
Аннотация
Работа посвящена решению важных задач подземной гидромеханики - задач многофазной фильтрации жидкостей при наличии скважин. Для решения этих задач представлены новые алгоритмы, основанные на методах декомпозиции. Алгоритмы реализованы на вычислительных системах нового поколения - гетерогенных суперкомпьютерах, построенных на основе современных центральных процессоров и графических ускорителей.
Ключевые слова: методы декомпозиции, многофазные течения, гетерогенные вычислительные системы.
Введение
Исследование процессов разработки месторождений углеводородного сырья с использованием математических моделей течений многофазной жидкости в пористых средах со скважинами является актуальной задачей. Математические модели таких процессов представляют собой системы связанных нелинейных нестационарных уравнений с частными производными. Для численного решения этих систем наиболее часто применяются вычислительные алгоритмы, основанные на тех или иных модификациях метода IMPES (неявная схема по давлению, явная - по насыщенности) [1]. Как правило, для решения задач многофазных течений в пористых средах используются сетки, сгущающиеся к интервалам вскрытия скважин [2, 3]. Требования к точности решения приводят к необходимости использования сеток, имеющих десятки миллионов узлов. Повышенные требования к вычислительным ресурсам обусловлены также нестационарностью и трехмерностью рассматриваемых процессов, неоднородностью физических характеристик пласта и др. Такие задачи могут быть решены только на высокопроизводительных вычислительных системах. Поэтому необходимы алгоритмы, позволяющие осуществлять распараллеливание. Наиболее популярными методами являются методы декомпозиции, когда производится разделение задачи на подзадачи, которые могут решаться параллельно. Одним из преимуществ методов декомпозиции является возможность их реализации на многопроцессорных вычислительных системах. Идея декомпозиции использовалась и ранее, но не в связи с параллельными вычислениями, и привела к появлению методов подструктур, подконструкций, макроэлементов, суперэлементов, фрагментов, модуль-элементов, редуцированных элементов, а также методов Шварца, матриц емкости и т. д. Такие методы всегда использовались как методы, позволяющие сводить решение исходной задачи в области со сложной границей к последовательности задач в подобластях, границы которых достаточно просты. С развитием гетерогенных вычислительных систем основной принцип разделения
71
72
А.В. ЦЕПАЕВ
области решения на подзадачи сводится к независимому решению этих подзадач на различных вычислительных комплексах (ядра, процессоры, графические ускорители). Это приводит к значительному уменьшению времени решения исходной задачи. При решении задач многофазной фильтрации на каждом временном шаге приходится определять поля давления и насыщенности. Для решения сеточных систем уравнений по давлению и насыщенности в работе [4] были предложены два различных метода декомпозиции области: один метод для решения сеточных уравнений по давлению, другой - по насыщенности. Метод декомпозиции области для определения поля давления основан на независимом решении систем алгебраических уравнений для сгущающихся участков сетки в подобластях и эффективном способе согласования этих решений с решением на грубой сетке. Для решения уравнения по насыщенности разработан новый метод декомпозиции области, основанный на сочетании элементов явной и неявной схем. На основе предложенных методов декомпозиции построены алгоритмы для решения задачи трехфазной фильтрации жидкости на гетерогенных вычислительных системах. Сочетание вычислительных ядер центрального процессора и графических устройств позволило значительно сократить время решения задач. Данные методы декомпозиции применялись также при решении задач фильтрации однофазной жидкости, подчиняющейся нелинейному закону Форхгеймера [5], напорно-безнапорной фильтрации [6].
1. Постановка задачи трехфазной фильтрации
Рассматривается трехфазная изотермическая фильтрация нефти, воды и газа, подчиняющаяся линейному закону Дарси. Считается, что пласт, нефть и вода несжимаемы и отсутствует массообмен между нефтяной и газовой фазой. Тогда справедлива следующая система уравнений:
qa
dSw
->+тчй
dS,
div(qw) + m div (q0) + m
o
dt
div (Pg qg) + :
d(Pg Sg)
g4g' ' dt
-(fak/pa)grad(p0
0,
0,
= 0,
- Paдz),
Y.Sa = 1,
a
pcao = pa po.
(1)
(2)
(3)
(4)
(5)
(6)
Здесь а соответствует одной из трех фаз: “о” (нефть), “w” (вода), “д” (газ); pa(x,y,z) - давление в фазе а, qa - вектор скорости фильтрации фазы а, Sa -насыщенность пласта фазой а, fa - относительная фазовая проницаемость, к -абсолютная проницаемость, ра - динамическая вязкость фазы а, pcao - капиллярное давление, Pa - плотности фазы а .
Область решения представляет из себя многосвязную область, внутренние поверхности которой определены поверхностями скважин в интервалах вскрытия пласта. В начальный момент времени считаются известными распределения давления и насыщенностей в пласте. На внешней поверхности пласта задаются граничные условия 1-го или 2-го рода. На скважинах задается либо забойное давление, либо суммарный расход жидкости при некотором фиксированном давлении, определяемом в процессе решения задачи жидкости. Таким образом, система уравнений
РЕШЕНИЕ ЗАДАЧ ТРЕХФАЗНОЙ ФИЛЬТРАЦИИ...
73
трехфазной фильтрации (1)—(6) дополняется граничными условиями
Pw = Рг на Г1, (7)
-(Ko + Kw + Kg) = qn на Г 2, (8)
Pw | dvk = Pk, к = 1, 2 - M, (9)
Sw Swr на Г3, So Sor на Г3? Sw Swk, k 1, 27 • • • , (10)
и начальными условиями
Pw = pw , Sw = sw, So = So0 в D, (11)
где riU Г2 = Г - внешняя граничная поверхность области D, Г3 - часть поверхности Г, через которую жидкость поступает в пласт, dVk - поверхность интервала вскрытия пласта к -й скважиной, Pk - заданное давление на к-й скважине, N -число скважин, M - число нагнетающих скважин Swk - заданная насыщенность в нагнетательной скважине.
2. Построение сетки
Системы дифференциальных уравнений аппроксимировались методом контрольных объемов. Построение сетки в методе контрольных объемов проводилось по следующему алгоритму.
1. Вся область решения покрывается грубой сеткой, элементы которой представляют собой параллелепипеды. Элементы построенной грубой сетки образуют множество E.
2. Определяется множество A из элементов множества E, имеющих общие точки с интервалами вскрытия скважин. Множество B состоит из элементов, примыкающих к границам множества A.
3. Элементы множеств A и B разбиваются вертикальными плоскостями по длине и ширине на 4 равных параллелепипеда. В случае, когда размеры полученных параллелепипедов по вертикали превосходят размеры по горизонтали, элементы дополнительно разбиваются горизонтальной плоскостью на 2 равных параллелепипеда.
4. Полученное множество параллелепипедов определяется как множество E , и переопределяются множества A и B .
5. В случае, когда размеры ячеек из множества A имеют размеры радиуса скважины, переходим к следующему пункту, в противном случае переходим к п. 3.
6. Из множества E исключаются элементы, отстоящие от поверхности интервалов вскрытия скважин на расстоянии менее радиуса скважин. Полученная область является многосвязной с внутренними границами, состоящими из граней элементов, окружающих интервалы вскрытия скважин.
7. Множество E дополняется криволинейными элементами, одна из граней которых является гранью внутренней границы. Противоположные грани определяются проектированием вершин граней внутренней границы на интервалы вскрытия по нормали к последним. На рис. 1 приведен криволинейный элемент в разрезе горизонтальной плоскостью.
Сетка Q определяется как объединение множества центральных точек элементов множества E и множества центральных точек граней внешней и внутренних поверхностей области, образуемой элементами множества E. Требуется определить сеточную функцию p = ра в узлах Q и значения насыщенностей Sa, в ячейках множества E. Расходы Qj, выходящие из j-й ячейки в i-ю, запишем в виде
74
А.В. ЦЕПАЕВ
интервал вскрытия ---—^ оси проекции
грань внутренней границы
Рис. 1. Криволинейный элемент в разрезе горизонтальной плоскостью
Qji = (Pi - Pj)/Rj , где Rij = Aij/(KWi + Koi + Kgi), Aj - коэффициенты, учитывающие сопротивление потоку жидкости между i-й и j-й ячейками за счет i -й ячейки, pi - давление в i-й ячейке. Для параллелепипедов имеем Aij = Lj/Dij, где Dij - площадь общей граничной поверхности i-й и j-й ячеек, Lij - расстояние от узлового значения i-й ячейки до общей граничной поверхности. Для криволинейных элементов коэффициенты Aij определялись численно с учетом движения потока по нормали к интервалу вскрытия. Сопротивление потоку жидкости в расходах через грани внешней и внутренних поверхностей учитывалось только за счет приграничных ячеек. Полученные расходы использовались для составления балансовых уравнений по давлению.
3. Решение задачи методами декомпозиции области
Задача решается в области D, представляющей собой пласт, ограниченный кровлей, подошвой, боковыми поверхностями и поверхностями интервалов вскры-
N
тия скважин Vk, к = 1, 2,..., N. Объединение |J Vk является дополнением мно-
k=1
госвязной области D до односвязной области. Область D разобъем на подобласти Dk удовлетворяющим условиям:
D = Do U (U Dk^ , Di П Dj = 0, i = j, Do П Dk = Yk, Sk П Yk = 0,
где Dk - прискважинные подобласти, Sk - суммарная поверхность интервалов вскрытия к-й скважины. На области D введем сетку Q, ячейки которой в прискважинных зонах Dk уменьшаются в размерах к интервалам вскрытия скважин по убывающей геометрической прогрессии. Для определения p^+1 на (n + 1)-м временном шаге используется алгоритм, подробно описанный в [4], с учетом капиллярных сил. Опишем кратко метод вычисления p^+1, основанный на методе декомпозиции области. Давления p^+1 вычисляются из системы (1)—(11), со значениями K0i, KWi, кп.
Для решения задачи (1)—(11), предлагается следующий метод декомпозиции. Вводятся дополнительные грубые сетки ид , покрывающие области Dk U Vk . Размер ячеек грубой сетки имеет тот же порядок, что и размер ячеек Q для Do. Решение p^+1 в области Do представляется как решение задачи фильтрации (1)-(6) P0+1 с краевыми условиями P0+1 = POt-yk на границах раздела Yk области Do и
прискважинных подобластей Dk , а в областях Dk - как сумма двух решений рО+к и p^k1. Решения p^k определяются на сетке Q в областях Dk, а решения pO+k -на сетках ^k в областях Dk U Vk, при этом р^^ = P^^ + pOa^ на границах Yk . Строится итерационный процесс, в котором на каждой i -й итерации независимо определяются решения рОак’1 в подобластях Dk, а затем совместно в области Do
на грубой сетке Q ив подобластях Dk на грубых сетках ид определяются ре-n+1’i
шения pa1k
в Dk
n+1’i n+1’i
и Pa3k . В построенном вычислительном процессе значения pa3k
стремятся к 0 и общее решение определяется рОкР
гл n+1’i
в Do и pa2k
РЕШЕНИЕ ЗАДАЧ ТРЕХФАЗНОЙ ФИЛЬТРАЦИИ...
75
Для определения насыщенностей Sa используется метод декомпозиции, описанный в [7], c учетом капиллярных сил, состоящий из следующих этапов.
1. Вычисляются полные расходы всех фаз, выходящие из ячеек грубой сетки в единицу времени
Qgr = (pj1 - pH 1 + Pag(Zj - Zi) /Rgj , (12)
где = Lij/Dij(KWi+K0i+Kgi)+Lii/Dij(KWj+K0j+Kgj), Dij - площадь общей граничной поверхности i-й и j-й ячеек, Lij - расстояние от узлового значения i -й ячейки до общей граничной поверхности, pgi - давление в фазе а в i-й ячейке.
2. Для полных расходов, выходящих из ячеек грубой сетки, вычисляются фазовые расходы по явной схеме (индекс jv)
Qgtljv = (Kg/(Kg + ко + Kg) up Qj,
где
up
(Kg/(KW + Kg + Kg))jt
Kgi/(KWi + Kgi + Kggi), Pai - Pg9Zi > Pgj - Pg9Zj,
Kaj/ (Kwj + Koj + Kgj ), pai — Pa9Zi < paj — Pa9Zj ■
3. Для каждой прискважинной зоны независимо вычисляются насыщенности по неявной схеме (индекс njv) из системы уравнений
mi
Vi (sgr1
Sgi)/At = £ Q
n+ 1,njv
aj,i
j
(13)
где сумма берется по j-м ячейкам, окружающим i-ю ячейку, Qgj1 = Qgji’jv для j -х ячеек грубой сетки, из которых расход поступает в прискважинную зону, Qgji ’njv = (Kg/(Kg + Kg + Kg))UpQgjl в остальных случаях. Фазовые расходы Qg+ijV являются граничными условиями при решении систем (13). Насыщенности
для ячеек грубой сетки Sgr1
g
Sgi
+ At/(miVi)Z2 Qj, где Qj
,n+1, njv
gj i
j
берутся из решения системы уравнений (13) для j -х ячеек прискважинных зон, из которых расход поступает в ячейки грубой сетки, QgГ = QgГ ’jv в остальных случаях. 4
4. Численные результаты
Предложенный алгоритм тестировался путем решения модельной трехмерной задачи трехфазной фильтрации жидкостей с различным числом вертикальных добывающих и нагнетающих скважин. Рассматривался десятислойный пласт (103 м х 103 м х18 м) с толщинами слоев d1 = 1 м, d2 = 1 м, d% = 3 м, d,4 = 1 м, d5 = 1 м, ер = 1 м, dg = 2 м, ds = 1 м, dg = 2 м, dm = 5 ми абсолютными проницаемостями k1 = 10-3 дарси, k-2 = 10-2 дарси, кз = 25 • 10-3 дарси, k4 = 10-3 дарси, k5 = 10-3 дарси, kg = 10-2 дарси, kg = 5 • 10-3 дарси, ks = = 10-2 дарси, kg = 10-3 дарси, k10 = 15 • 10-3 дарси соответственно. Кровля пласта считалась непроницаемой, на боковых поверхностях и подошве пласта давление рг = 12 МПа, на скважинах Pk = 3 МПа, на боковой поверхности насыщенность Sw =0, на подошве Sw = 1 .В начальный момент времени So = 1, Sw =0, за исключением расположенной в центре пласта у кровли зоны Гд, занимаемой газом, в которой So = 0, Sw = 0, Sg = 1. Динамическая вязкость воды ^w = 1 мПа • с, динамическая вязкость нефти ^o = 15 мПа • с, плотность нефти po = 882 кг/м3,
76
А.В. ЦЕПАЕВ
плотность воды pw = 103 кг/м3. Относительные фазовые проницаемости полагались линейными функциями от насыщенностей, функции pcao - логарифмическими функциями от насыщенностей [1], Pcow = Pcow (Sw) , Pcow = Pcow (Sw), pcog = pcog (Sg). Каждый интервал вскрытия моделировался круговым цилиндром с радиусом основания r = 0.1 м и замыкался сверху и снизу сферическими поверхностями радиуса r = 0.1 м. Таким образом, для каждой точки поверхности интервалов вскрытия вектор нормали определен однозначно. Ячейки, примыкающие к скважинам, имели размеры как в горизонтальной плоскости, так и по высоте, порядка 0.1 м.
На основе предложенных методов построены алгоритмы для решения задачи на суперкомпьютерах с большим числом вычислительных ядер и графических процессоров. Алгоритмы распараллелены с помощью MPI-процессов, OpenMP- и CUDA-технологий. Использовались язык С+—Н и среда разработки приложений Visual Studio 2008. При решении задачи для каждого MPI-процесса выделялось равное число сгущающихся участков сетки. Для решения задач на сгущающихся участках, соответствующих одному MPI-процессу, порождались потоки с помощью технологии OpenMP, которые распределяли эти задачи на ядра процессора и на графические устройства. При этом задачи распределялись динамически, то есть по мере их решения. Вычисления на графическом устройстве включали в себя:
1) построение систем уравнений для определения поля давления и насыщенности;
2) решение линейных систем уравнений для определения поля давления с помощью библиотек Cublas и Cusparse;
3) решение нелинейных систем уравнений для определения полей насыщенностей в прискважинной зоне;
4) определение полей насыщенности по явным схемам для ячеек внескважинной зоны.
Предложенные алгоритмы тестировались на кластере, состоящем из 4-ядерных вычислительных узлов с процессорами Intel Core i7 2600, с графическими ускорителями компании NVidia GTX 560 TI путем решения модельной трехмерной задачи трехфазной фильтрации жидкости с вертикальными добывающими и нагнетательными скважинами. Каждая прискважинная зона содержала 25000 узлов, внескважинная зона - 6000 узлов. Общее число узлов для 200 сгущающихся участков достигало 5 х 106. Задача решалась на двух вычислительных узлах с общим числом ядер 8 (4 ядра на узел) и 2 GPU-ускорителя. Задачи запускались в режиме 1 MPI-процесс на 1 узел, каждому MPI-процессу соответствовало 1 GPU-устройство. Рассматривались следующие варианты запуска задач.
1. На всех доступных ядрах с использованием одного MPI-процесса и без использования GPU-ускорителей.
2. На ядрах с использованием MPI-процессов и GPU-устройств, равных числу компьютеров в кластере. В этом случае каждому MPI-процессу соответствовало одно GPU-устройство и одно ядро.
3. На всех доступных ядрах с использованием MPI-процессов и GPU-устройств. В этом случае использовались все доступные ресурсы кластера (2 GPU-устройства и 8 ядер).
В табл. 1-3 приведены результаты решения. Сравнение проводилось с решением задачи на одном ядре без использования GPU-устройств.
В табл. 1 показано ускорение времени решения задачи, при использовании ядер только одного вычислительного узла. Эффективность использования ядер достигала 72%. Эффективность использования GPU-устройств составила 80% (см. табл. 2). Предложенными методами решалась также задача, в которой область,
РЕШЕНИЕ ЗАДАЧ ТРЕХФАЗНОЙ ФИЛЬТРАЦИИ...
77
Табл. 1
Вариант 1
Число ядер 1 2 3 4
Ускорение 1 1.7 2.1 2.9
Табл. 2
Вариант 2
Число GPU-устройств 1 2
Ускорение 17.3 30.1
Табл. 3
Вариант 3
Число GPU-устройств 1 2
Ускорение 20.2 36.3
занимаемая газом, предполагалась пустой, то есть фактически рассматривалась двухфазная фильтрация [8]. Задача решалась на суперкомпьютере «ГрафИТ!» Научно-исследовательского вычислительного центра Московского государственного университета ГУ имени М.В. Ломоносова. На момент запуска задачи на суперкомпьютере «ГрафИТ!» было доступно 5 вычислительных улов с общим числом ядер 60 (12 ядер на узел) и 15 GPU-ускорителей. Так же, как и в предыдущем примере, рассматривались 3 варианта запуска задач. На рис. 1-3 приведены результаты решения.
На рис. 2 проиллюстрировано ускорение времени решения задачи в зависимости от числа ядер в случае запуска задачи на одном узле с одним MPI-процессом и без использования GPU-устройств (1 вариант). Сравнение проводилось с решением на одном ядре и без использования GPU-устройств. Из рисунка видно, что на 12 ядрах получено семикратное ускорение (эффективность 60%). На рис. 3 приведено ускорение времени решения задачи в зависимости от числа MPI-процессов в случае запуска задачи с использованием различного числа GPU-устройств и выделением одного ядра для каждого MPI-процесса (2 вариант). Сравнение проводилось с результатом решения на одном ядре и одном GPU-устройстве. Достигнута эффективность 70% при решении на 15 GPU-процессоров.
На рис. 4 приведено ускорение времени решения в случае запуска задач с различным числом MPI-процессов и GPU-устройств (вариант 3), при этом использовались все доступные ядра узлов, на которых запускались MPI-процессы. Сравнение проводилось со временем решения задачи на одном GPU-устройстве и одном ядре. Отметим, что если сравнивать решение на одном ядре без MPI-процессов и GPU-устройств с решением задачи с использованием всех доступных вычислительных устройств суперкомпьютера, то ускорение достигало 120 раз.
Заключение
Построены алгоритмы для решения задачи трехфазной фильтрации жидкости на сетках со сгущениями в районе скважин с применением методов декомпозиции области. Декомпозиция сеточной системы уравнений по давлению основана на согласовании решений для сгущающихся участков с решением на грубой сетке за счет введения дополнительных грубых сеток. Декомпозиция сеточной системы уравнений по насыщенности осуществляется путем независимого решения уравнений на сгущающихся участках по неявным схемам и согласования этих решений с решением на грубой сетке с использованием элементов явной и неявной
78
А.В. ЦЕПАЕВ
1 2 3 4 5 6 7 & 9 10 11 12
число ядер
Рис. 2. Ускорение времени решения задачи на одном узле в зависимости от числа ядер (вариант 1)
123456789 101112 131415
число MPI процессов
Рис. 3. Ускорение времени решения задачи в зависимости от числа MPI-процессов (вариант 2)
Рис. 4. Ускорение времени решения задачи в зависимости от числа MPI-процессов (вариант 3)
схем. На основе методов декомпозиции созданы алгоритмы для решения задачи на гетерогенных вычислительных системах. Установлена высокая эффективность использования многопроцессорных систем, построенных на базе графических процессоров.
Работа выполнена при поддержке РФФИ (проект № 12-01-31180).
РЕШЕНИЕ ЗАДАЧ ТРЕХФАЗНОЙ ФИЛЬТРАЦИИ...
79
Summary
A. V. Tsepaev. Solving of Three-Phase Flow Equations with Hydrodynamically Imperfect Wells on Heterogeneous Computing Systems.
The work is aimed to solve the important problems of underground fluid mechanics, i.e. the multiphase fluid filtration problems. To solve these problems, new numerical algorithms based on domain decomposition methods have been developed. These algorithms have been realized on new generation heterogeneous computing systems, built using modern central processing units and graphics accelerators.
Keywords: decomposition methods, multiphase flows, heterogeneous computing systems.
Литература
1. Chen Z., Huan G., Ma Y. Computational methods for multiphase flows in porous media. -SIAM, 2006. - 557 p.
2. Басниев К.С., Власов А.М., Кочина И.М., Максимов В.М. Подземная гидравлика. -М.: Недра, 1986. - 303 с.
3. Азис Х., Сеттари Э. Математическое моделирование пластовых систем. - М.: Недра, 1982. - 407 с.
4. Мазуров П.А., Цепаев А.В. Алгоритмы для распараллеливания решения задач двухфазной фильтрации жидкости на сетках со сгущающимися участками // Вычисл. методы и программирование. - 2006. - Т. 7, № 2. - С. 115-123.
5. Мазуров П.А., Цепаев А.В. Метод решения нелинейных задач фильтрации жидкости в трехмерных пластах с гидродинамически несовершенными скважинами // Матем. моделирование. - 2004. - Т. 16, № 3. - С. 33-42.
6. Губайдуллин Д.А., Мазуров П.А., Цепаев А.В. Алгоритм решения трехмерных задач напорно-безнапорной стационарной фильтрации жидкости со сгущающимися участками сетки // Вычисл. методы и программирование. - 2005. - Т. 6, № 2. - С. 217-225.
7. Губайдуллин Д.А., Никифоров А.И., Цепаев А.В. Алгоритмы распараллеливания на сгущающихся сетках в задачах трехфазной фильтрации жидкости // Вычисл. методы и программирование. - 2007. - Т. 8, № 2. - С. 360-366.
8. Цепаев А.В. Использование гетерогенных вычислительных систем для решения задач фильтрации жидкости методами декомпозиции области // Вычисл. методы и программирование. - 2012. - Т. 13, № 1. - С. 39-44.
Поступила в редакцию 28.06.13
Цепаев Алексей Викторович - кандидат физико-математических наук, старший научный сотрудник, Институт механики и машиностроения КазНЦ РАН, г. Казань, Россия.
E-mail: alex_ tsepaev @mail. ru
