CC BY
15УДК 519.6/.87:532:622.276:004.9
ПОЛНОСТЬЮ НЕЯВНАЯ СХЕМА РЕШЕНИЯ ЗАДАЧ ТРЕХФАЗНОЙ ФИЛЬТРАЦИИ НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ В ПАКЕТЕ ПРОГРАММ НИМФА
О. И. Бутнев, И. В. Горев, С. С. Колесников, В. Ю. Кузнецов, В. А. Пронин, М. Л. Сидоров,
А. Д. Яруллин
Российский федеральный ядерный центр — Всероссийский научно-исследовательский институт
экспериментальной физики, boi@vniief.ru
В статье представлен обзор отечественного программного комплекса НИМФА. Данный программный комплекс позволяет проводить расчеты задач многофазной фильтрации по модели «Black Oil» на высокопроизводительных системах с распределенной памятью. Представлены особенности ПК НИМФА: полностью неявная схема решения уравнений модели «Black Oil» на неструктурированных сетках, параллельный сеточный генератор, принципы распараллеливания вычислительных модулей. Для данного ПК проведены различные тесты из набора SPE. Получено хорошее качественное и количественное согласие результатов с коммерческими программными продуктами.
Ключевые слова: многофазная фильтрация, высокопроизводительные системы, неструктурированные сетки, модель «чёрной нефти», параллельный сеточный генератор, полностью неявная схема.
A FULLY IMPLICIT SOLUTION PATTERN FOR THREE-PHASE FILTERING PROBLEMS APPLIED TO THE UNSTRUCTURED GRIDS IN THE NIMFA SOFTWARE PACKAGE
O. I. Butnev, I. V. Gorev, S. S. Kolesnikov, V. Yu. Kuznetsov, V. A. Pronin, M. L. Sidorov,
A. D. Yarullin
Russian Federal Nuclear Center — All-Russian Research Institute of Experimental Physics, boi@vniief.ru
The article presents a review of Russian application software NIMFA. This application software allows calculations of problems of multiphase filtration model "Black Oil" high-performance systems with distributed-memory. Features of application software NIMFA are fully implicit scheme the solutions of the equations of the model "Black Oil" on unstructured meshes, parallel mesh generator, the principles of parallel computing modules. For this application software carried out the various tests from the set of SPE. Obtained good qualitative and quantitative agreement of the results with commercial software products.
Keywords: multiphase flows, high performance, unstructured mesh, model "Black Oil a parallel mesh generator, fully implicit scheme.
Введение
Разработка месторождений углеводородов представляет собой комплексную проблему, для успешного решения которой требуется привлечение знаний и опыта, накопленных в различных областях науки и инженерной практики. Объекты нефтегазовой отрасли имеют ряд особенностей: большие размеры по площади, достигающие сотен квадратных километров, сложную структуру пластов и сложные физические процессы, проходящие при их разработке. Также при решении задач данной отрасли дополнительную трудность вызывает разномасштабность моделируемых объектов. Приведенные выше особенности накладывают специфические требования к математическим методам их моделирования.
На данный момент для моделирования процессов фильтрации в геологических пластах применяются в основном зарубежные коммерческие программные продукты, такие как ECLIPSE компании Schlumberger [1], Tempest More компании Roxаr [2], STARS компании CMG [3]. Российские программные продукты в практике математического моделирования процессов нефтедобычи практически не используются. Зарубежные программные продукты имеют, во-первых, дорогостоящую лицензию и, во-вторых, в условиях экономических санкций, введенных против Российской Федерации, нет гарантий, что в ближайшей перспективе не будет введено эмбарго на поставку оборудования и объектов интеллектуальной собственности иностранного происхождения. В связи с этим остро встает вопрос о разработке отечественных пакетов программ для решения задач фильтрации с применением технологий суперкомпьютерного моделирования.
Необходимость внедрения современных методов математического моделирования определяется существенным ухудшением структуры запасов нефти и газа в РФ. Вследствие этого возникает задача, связанная с повышением нефтеотдачи, для решения которой необходима разработка принципиально новых технологий воздействия на пласт при добыче нефти, которые могут быть смоделированы только с использованием новых быстродействующих суперкомпьютеров. Следует отметить, что отработка данных технологий путем только практических экспериментов является весьма дорогостоящей, а получаемые результаты при таких тестах очень ограничены.
Суперкомпьютерные технологии - один из ключевых факторов конкурентоспособности отечественной нефтегазовой отрасли. Для обеспечения устойчивого экономического развития России требуется создать технологический паритет отечественных предприятий нефтегазовой отрасли с лидерами мирового рынка. Важной сферой применения математического моделирования задач фильтрации является решение проблем прогнозирования, контроля и управления процессами разработки пластов с целью повышения их нефтеотдачи, - в этом состоит и основное коммерческое использование программных продуктов. Одним из таких программных комплексов является комплекс НИМФА [4,5], разрабатываемый во ФГУП «РФЯЦ-ВНИИЭФ».
На данный момент в ПК НИМФА реализована модель трехфазной фильтрации «Black Oil» (модель «черной нефти»). Данная модель покрывает существенный класс двух и трехфазных задач [6]. Для дискретизации по пространству используется метод конечных объемов. Для аппроксимации потоковых слагаемых используется метод «отложенной» коррекции [7].
Математическое описание модели «Black Oil»
Модель «Black Oil» (модель нелетучей нефти) содержит три фазы: нефть (о), воду (w) и газ (g). Вода и нефть не смешиваются между собой и не обмениваются массами. Газ предполагается растворимым в воде и нефти [15].
Уравнения нелетучей нефти представляют собой систему уравнений, объединяющей уравнения сохранения массы и уравнения движения Дарси [6],[8]. Движение среды описывается с точки зрения Эйлера. Система дифференциальных уравнений сохранения массы, описывающая трехфазное трехмерное течение жидкости в пористой среде, имеет вид [15]:
д_ dt
т
д_ dt
{ (f>Sw\ V Bw )
+ V
+ V
о
BO \ /
W
w
Bw
\
— qo;,
— qw ;
9 ((t)Sg (f>RsoSo +
д V B.
g
Bo
<f>RswSw Bw
+ V
—»•
W RsoW g- + +
R
—»•'
sw W w
\
B
g
Bo
Bw
— qoRso + qwRsw + qg •
(1)
(2)
(3)
/
Уравнения движения (скорость фильтрации) описывается законом Дарси:
Wo — - Ао (Vpo- 7oVh), —»•
Ww — - Aw (V pw- 7wVh), —»•
Wg — Ag(Vpg-7gVÄ) •
Здесь Ao, Aw, Ag — коэффициенты проводимости соответствующей фазы — имеют вид:
, _ kmKa ^ krwKa ^ krgKa Ao — ; Л w — ; ^g — •
(4)
(5)
(6)
ßo
w
g
Здесь ф — пористость породы [—]; Ка — тензор абсолютной проницаемости породы [Ь 2]; £Г0, кгё — относительные фазовые проницаемости нефти, воды и газа[—]; //0, цш, % — коэффициенты
k
динамической вязкости нефти, воды и газа [МТ~ ']; ро, рт, рё — давление, соответственно, в нефтяной, водяной и газовой фазах [МТ~2Ь_']; 7о = роё,^(т = ртё, 7ё = Рёё — веса соответствующих фаз, ё — ускорение свободного падения [ЬТ~2]; Н — глубина относительно уровня моря, отсчет ведется вниз по вертикали [Ь], да — объемный источник или сток а-той фазы (а = о,т,ё).
Для замыкания системы уравнений (1)-(6) необходимы следующие соотношения:
ро ~ рт = Рсот (3о,3т) , (7)
рё - ро = Рсоё ($о,$ё) , (8)
So + 3т + = 1. (9)
Отношения (7)—(8) выражают зависимости капиллярных давлений в системе вода-нефть и газ-нефть, используются для пересчета давления воды рт и давления газар^ соответственно. Уравнение (9) используется для пересчета насыщенностей.
Для полной формулировки математической модели исходная система уравнений (1)—(9) должна быть дополнена уравнениями состояния (данными PVT), а также начальными и краевыми условиями.
Под начальными условиями понимается начальное распределение искомых значений давлений и насыщенностей на момент времени t = 0.
Под краевыми условиями понимается задание граничных условий (режимов работы) на границе моделируемой области и на каждом из источников, представляющих скважины.
В системе уравнений (1)—(9) 17 неизвестных:
ро зо Бо о кго Рзо
рт зт Бт Рт кгт Рзт
рё Зё Бё Рё кгё,
для которых требуется 17 уравнений:
- 3 уравнения состояния или таблицы (PVT), позволяющие вычислить ро,рт,рё = /(Р);
- 3 уравнения или таблицы PVT, задающие Бо,Бт,Бё = /(Р);
- 2 уравнения или таблицы PVT, задающие Рзо,Рзт = !(Р);
- 3 таблицы относительных фазовых проницаемостей, задающих кго = /(зо), кгт = /(зт), кгё = /(зё);
- 1 уравнение для насыщенностей, дающее зё = 1 — зт — зо;
- 2 таблицы капиллярных давлений, дающие рт = ро — Рсот ,рё = ро + Рсоё;
Таким образом, для нахождения оставшихся переменных: ро,зт,зо требуется 3 уравнения (1)-(3).
Метод учета скважин в численных моделях фильтрации основан на допущении того, что вблизи скважины течение описывается аналитическим решением, граничные условия для которого определяются из численного решения задачи для пласта. Этот подход впервые предложен Д. Писманом [9].
Скважина рассматривается как множество интервалов перфорации, вскрывающих сеточные блоки (ячейки). На скважине могут задаваться забойное или устьевое давление, дебит флюидов в поверхностных или пластовых условиях. Распределение дебита по фазам и по интервалам перфорации осуществляется с учетом гравитационных сил и состава смеси на скважине.
В случае, когда на скважине заданы забойные давления, объемный фазовый дебит интервала перфорации, приведенный к стандартным условиям, имеет вид [15]:
д1о = {КаН1 (В^) [р1-р*,ъ-1 {Н^ Нге{)] |; I = 1, 2,.., Мр;а = о,т,ё; (10) %Т = % + Х] ШаЪа;1 = 1, 2,..., Ыр;а = о,т; (11)
о
где р1 - давление в ячейке, содержащей интервал перфорации I; ртъ - забойное давление, приведенное к опорной глубине; - коэффициент скважины для интервала перфорации I; Н - толщина интервала I (может быть меньше, чем толщина ячейки, содержащей интервал I, но не может быть больше ее);
- индекс противотока; 7 - осредненный удельный вес жидкости в стволе скважины; Н[ - глубина интервала I; Нге{ - опорная глубина; Ыр - число интервалов перфорации; т - полный дебит газа.
Здесь для добывающей скважины должно выполняться соотношение:
6/ =
и
1, Р/^РшЬ^1 {Н/^ Н^) > О, о, Р/ ~ РшЬ ~ 7 (Н/- Нге{) < О,
в котором для нагнетательной скважины знаки неравенств меняются на противоположные; 7 =
Е/ Ео Ч/сЛ/а! Е/ Ео Ч/а, причем Е / Ео Ч/а = 0, в противном слу^е 7 = 7^ где 7а согласуется
с начальными условиями.
Коэффициент скважины в выражении (10) вычисляется как
ш/ = 2'
(г)+^
где должно выполняться условие Гь/ > гш/ • ехр (—3/); Гь/ - эквивалентный радиус сеточного блока, содержащего интервал V, Гш/ - радиус скважины на интервале V, 3/ - скин-фактор на интервале Ь
Эквивалентный радиус сеточного блока, содержащего интервал 4 рассчитывается с помощью формулы Писмана:
£_
Гь/ = 0,28—-=-—-,
4+ 4 Кл
V кх/ у ку/
где кх/ и ку1 - абсолютные проницаемости сеточного блока, содержащего интервал /, по соответствующим направлениям; Д X/ и А у/ - размеры сеточного блока, содержащего интервал I Фазовые дебиты скважин рассчитываются с помощью выражения
Чо = Ч/<* ■ (12)
/
В случае, когда на скважине заданы фазовые дебиты Ча, « = 0,ш^, либо дебиты жидкости Чь = Ч , = 0,ш, используется упрощенный метод распределения дебитов по интервалам перфорации:
В/аТ/аЧа
Ч^ = ^ р т ;« = 0,ш,я; (13)
/__</ В/а' /а
В/аТ/0Чь а ,лл.
Ч/а = ^ ^ р ~ ;« = 0,ш;Р = 0,ш,я; (14)
Ъ/ЪрЩ'/р
где 7/а = Ш1/к/Н/В^Кехр, аналогично для Т/р.
Для случая, когда на скважине заданы фазовые дебиты, либо дебиты жидкости, забойные давления рассчитываются из выражения (10).
Для случая, когда на скважине задан дебит жидкости, по формуле (14) производится распределение дебитов нефти и воды по интервалам перфорации. Затем по формуле (10), записанной для нефти или воды, определяется забойное давление. После чего по формуле (11) определяется дебит газа [15].
Искомые параметры скважины (забойное давление либо дебит) могут быть учтены в явной или неявной форме.
Сеточный генератор
Особый интерес представляют расчеты крупных нефтяных месторождений. При расчете больших территорий покрываемая площадь может составлять от сотен до нескольких тысяч квадратных километров. Например, при расчете месторождения площадью в 5000 кв.км. и 100 м. в глубину со средним размером ячейки (20м. х 20м. х 1м.) потребуется порядка 1.25 млрд. ячеек. Для расчета таких задач необходим высокопараллельный комплекс с параллельным генератором сеток.
Исследуемый объект представляет собой математическую модель коллектора, месторождения, участка почвы и т.п. Данная модель образуется исходя из геологической модели. Как правило, математические модели исследуемых объектов представляют собой слоистые модели, полученные экструзией.
Рис. 1. Пример сетки в клиновидном пласте Рис. 2. Пример сетки в клиновидном пласте
Рис. 3. Пример дискретизации русла реки Рис. 4. Фрагмент сетки
В пакете программ НИМФА используются неструктурированные сетки, удобные для описания сложных структурных элементов геологической модели с целью повышения точности расчета вблизи геометрических особенностей.
Авторами пакета НИМФА был разработан собственный параллельный генератор трехмерных неструктурированных пластовых сеток, покрывающий и превосходящий возможности зарубежных коммерческих аналогов как по качеству получаемой сетки, так и по скорости построения. Наиболее подробно сеточный генератор пакета НИМФА представлен в [5]. В генераторе используется интерфейс MPI. Основная суть метода построения сетки заключается в следующем:
1) строится планарная неструктурированная сетка в параллельном режиме,
2) на основе построенной планарной сетки строится объемная сетка путем экструзии с применением метода обобщенной угловой точки,
3) производится адаптация сетки к гидрогеологическим объектам.
С помощью разработанного сеточного генератора на примере одной из задач фильтрации жидкости в пласте со скважинами в параллельном режиме была получена трехмерная многогранная сетка в ячеечно-граневом представлении с числом ячеек более 1.1 миллиарда с использованием 21600 процессорных ядер. Время генерации сетки составило около 4 минут.
На рисунках 1-6 приведены примеры построения сеток методом обобщенной угловой точки.
Численная модель
В пакете НИМФА для получения разностных уравнений для модели «черной нефти» используется метод конечных объемов. Для аппроксимации потоков через боковые грани ячеек используется метод отложенной коррекции.
Для решения разностных уравнений используется полностью неявный метод, известный ещё как метод совместного решения [6, 8].
Суть данного метода заключается в записи производных от насыщенностей в правой части уравнений (1)-(3) через производные от давлений и решение уравнений в переменных давления нефти po, насыщенности нефти So и насыщенности воды Sw. Впервые это было предложено Дугласом [10]. Позднее Коутсом [11] было предложено использовать полностью неявный метод для моделирования многофазных течений в трещиновато-пористых средах.
Э г
а ■<>'■ ■
[ Г ':-,.- .' :
.........../-¡Л;' '..'":'.
тшм
а
Рис. 5. Пример сетки с дроблением во всех пластах
тт
шЯж
Рис. 6. Пример сетки с дроблением в некоторых пластах
Так как исходные уравнения являются нелинейными, в качестве метода линеаризации используется метод Ньютона—Рафсона [8].
В результате дискретизации и линеаризации получается по три уравнения для каждого сеточного блока. При этом матрица становится блочного вида. Для наглядности приведем пример. На рисунке 7 представлена плоская двумерная область. Структура блочной матрицы для такой области представлена на рисунке 8.
Ненулевые элементы помечены «х». Для решения данной матрицы будет использоваться библиотека параллельных решателей LParSol версии 2.0 [12].
Методы получения сеточных уравнений описаны ниже. Метод Ньютона—Рафсона
Пусть мы имеем общую систему нелинейных дифференциальных уравнений:
СтШр(х)]} = Ш), т = 1,2,....М, хе О, (15)
где £т обозначается линейный дифференциальный оператор, Бт(•) является нелинейной функцией,
Р = (р\ ,Р2,____,Рм)Т являются зависимыми переменными, / = (/ь/2 ,....,/м)Т — заданный вектор, и M
общее количество уравнений. С помощью метода Ньютона—Рафсона для уравнений (10) мы получим систему итерационных уравнений. Применяя разложение в ряд Тейлора для Рт(р + 8р), получим
Рт(р + 5р) = Рт(р) + УРт(р) • 5р + 0(|5р|2), (16)
где |<5р|является евклидовой нормой для 8р. Если слагаемыми высокого порядка 0(|5р|2)пренебречь, то Рт(р + 8р) может быть аппроксимировано как
Рт(р + Ьр)^ Рт(р) + УРт(р)> ¿р. (17)
Если подставить (12) в (10), мы получим итерационное уравнение
а?т(р1) + УРт(р1) • У] = Ш), т = \,....м, хе О, (18)
где р1 есть решение на 1-ой итерации вектора р и^Рт(р1) является УРт(р) при р = р1, с начальным решением р0. В итерационной системе (13) поправочный вектор 8р1 является неизвестным. Система может быть переписана в виде
)• 5р1 ] = Вт(х), т = 1,....м, х е О, (19)
где дт(х) = ¡т — £,[Рт(р1)], а\!¥т (р1 )и Рт(р1 )считаются фиксированными.
1 g 11
1 5 з 4
5 6 7 8
9 10 11 12
Рис. 7. Двумерная сеточная область
Рис. 8. Структура матрицы для сеточной области, показанной на рисунке 7, при использовании неявного метода
Новый вектор решения р1+1 получается добавлением поправочного вектора к решению с предыдущей итерации, то есть
р1+\ = р1 + $р1.
Этот итерационный процесс сходится при уменьшении 6р1 ниже заданного е. Метод совместного решения
Пусть п > 0 обозначает временной шаг. Введем следующий временной оператор для функции от времени V:
8у = vn+\ - V'1'.
При неявной аппроксимации по времени система уравнений (1)—(3) определяется как
5
(BW) -с1т + qw+1,
^ (BO) fo+^h)]+ qo,
+ nBr) =A^+1(AP°+1- h)] +
+ Д[Аno+lR0 + '(ApO+1 - lno+X&h)] + qg+1 + RO+lq0+1,
(20)
где At = tOtO. Система (15) является нелинейной с неизвестными pO+1 и SO+1, a = o,w,g, и может быть линеаризована методом Ньютона—Рафсона. Для этого мы запишем
POO+l,l+1 = POO+Ul + p, SO +1J+1 = SO +u + S, a = o,w,g,
где l обозначает номер итерации, а Spa и 5Sa представляют изменение давления и насыщенности за этот итерационный шаг. Зная, что
Vo+1 ^ vO+1,1+1 = уO+l,l + gy
Sv И vO+1,1 -VO + Sv, перепишем систему (15) в следующем виде:
1
Ы \
ы \
ы
\Вп) \Вп ) \Вп )
(£).....Ч^Г
= д[л П+^+Чд рп+\'м- Н)] + ч1+\'м
= А[Хп0+\'/+\{Ар'+\'1+\ - Ц+\'1+\ДН)\ + ц'п+\'/+\
\Вё Во ) Х^Вё Во ) \вё Во )
= ^[>!ёп+\'/+\(^рПп+\'/+^ 1пё+\'/+\АН)]+А[\п0+\'/+\рп+\'/+\(Арп+\'/+\^ -0+\-/+^Н)] + + +\./+\+ ^+\'/+\д'+\'/+\,
(21)
В этой системе 5р0 и являются неизвестными. В нашем случае основными переменными являются
8р, 6 Б
Отметим, что = — 6Бо — 6Бп. Соответственно, левая часть системы (21) может быть представлена следующим образом. Для водной фазы
6 (ФБпЛ
\Вп)
— Спр$р + СпЗп$БП;
где
Для нефтяной фазы
спр = Ф%(ВВп) + (^-^г) , ^ = (£) .
5 (^В^) = Сор^р + СоБ^Бо;
где
Сор — МШ +(«Т) СоБо = Ц)
Для газовой фазы
ц Б'+ъл Вё во
(
— С'рдр + С'Зп$Бп + С'Б^Бо;
где
сёр =
10^ + ВзоБЛ/ + ■ Ве Во )
/
-В'
Ф1Бё-Вг +
и;
С'Бп ^ С'Бо ^ В°).
(22)
(23)
(24)
Опуская далее временной индекс, запишем значение проводимости фазы на / + \ итерации
х/+\ = ° р/+\В/0+\
К, а = о,п,ё,
где рП+\ = рп. Таким образом, дебиты скважин определяются следующим образом
Мп мп
дП+\ = Е Е ™(ит) ^ Ы2)ж - р/+\ + РС+П- р1п+хр(гь1~ *)) *(*-
. , , Рш \ 7
¡=\ т=\
м„ Щ Ь1+х
= Е Е ШМ ~Т+1 {(Рьн)М- рМ- Р1о+1Р(*Ь1- *)) 6(х- х(1т)), (26)
1=1 т=1 V0
ми+1
?Г=ЕЕ Ш (1т) % р1+{ + р+ - ^Р^ г)) Кх- х(М).
1=1 т=1 %
Для удобства введем понятие фазового потенциала
ф 1+1 = рI рI + 1 н
■^т С Ссот !т у "•>
Ф0+1 = р1н,
фМ= р1+1 + рС+1- Н.
Далее мы представим потенциалы, проводимости и дебиты фаз в основных переменных: о и 8Бт.
Для водной фазы
Ф1т+1 = Ф т + йтрбр + йтБт^т, (27)
где
—тр — 1
[—т)^ ^ = {
Для нефтяной фазы
Ф0+1 = Ф о + ¿орбр, (28)
где
&ор = 1 " Ш рН.
Для газовой фазы
Л^1 — Л1 _1_ И Хп _1_ И „(XI + ,, ,
)т т ооо
Ф^1 =ф1е + йёр8р + (8Бт + 6Бо), (29)
где
^р = ^ ( ) рН &орь =~ рН
где
1 (^соЛ1
V )
Переносимости фаз могут быть выражены похожим образом. Для водяной фазы
тт+ = тт + Етр8р + ЕтЗт$Бт, (30)
Е =( Кт—Б^Л1 к Е =( й^ш 1 \1 к
Етр — I , I к, ЕтБт — I ,с 0 I к.
\11т йр } \ йБт ¡ЛтВт}
Для нефтяной фазы
то+1 = то + Еордр + ЕоБт5Бт + ЕоБ^Бо, (31)
где
Ео = (кго—(—1V К Ео5о =(——11 к Еш5ш =( (й^о^ 1 К
р \ йр \^оБо)) ' цоБ
о / \ \ —Бт Б
Для газовой фазы
т^1 = т1ё + Еёр8р + Еёз (8Бт + 8Бо), (32)
где
Е- = ^ = - ;) к.
Дебиты скважин также представляются подобным образом. Для водяной фазы
Мп Мп
яП+\ = чП + Е Е шЫт) [еПр$р + ПпБ + ПпбрЦ 8{х- хЫт)),
}=\ т=\
(33)
где
еп) =_±
епр = Мп
е(1) =
пБп Цп
\ + -р Р{гЬН
Р — -Тп
,епрЬН = Цп ,
—п {рЫ - р ~ рСоп - рпр{гЫ
-ра
1гп -Бп
Для нефтяной фазы
Мп Мп
Чо+\ = Чо + Е Е Ш{1,т) еЫ)*р + еоЦёБп + е^о + е^МП 8{х - х{Ы,т)), (34)
Ы=\ т=\
где
еЫ = { к
сор —
-Ио
-р
рЫ^ р-
рф{НЫ
еЫ) =
еоБп ~
е( )= еоБо =
Для газовой фазы
X [ -кго ¡Ло \ -Бп
-кго \ -Бё Цо
-кгс -Бд
^р^ РФ {НЬН
Ы
ЬН
Н))~\ П + --оР{Н
Ы)
М)
ЬН
Ыр-рорН0] еорьн = (£) .
Мп Мщ
чё+\= чё + Ыт)
Ы=\ т=\
еЫ>&р + + 5Бо) + е§рш8рЫ 6{х - х{Ыт))
(35)
где
еЫ = кгё
{к" (■
^р^ рСоё- РФ {НъН - Н
1 + -рр {НЫ
ё
еЫ) =
еёБ =
{рЫ)~ р- рСое- РёРЫ Н)) - кгё--Щ-)_
еорьн = I ■•
ё
Представим так же подобные разложения для коэффициентов растворимости Я/о и объемного расширения В0,а = о,ё,п.
\ = Х/о + Ърбр, В/+\ = В/{^ Ьар8р), а = о,п,ё,
]/+\ - о/
(36)
где
г/р —
Ш, Ьр = - "В)
р В р
а = о, п, ё.
В -р
Итоговые разностные уравнения
Подставив выражения (22)-(36) в уравнения (20), а также пренебрегая слагаемыми высокого порядка (типа 8р • 8Бо), мы получим систему линейных уравнений:
{В) (Вп) + Сп^р+
= У- {{тп + Епр8р + Бпз^Бп)
+ В\г{ч1п + VI{Ыт)
у#п) + V • {тпу-пр8р) + V • ТУ {-пБ^Бп)) + Ыр + е^Бп + е^ьнШ 6{х- х{Ы^)] + Ьпрч
для водяной фазы;
р
(37)
{Во)1 ^{В)" + Сор8р + СоЗо8Б>
= V • {{то + Бор8р + БоБп$Бп + ЕбоБ) VФ0) ^ • {Т/У {-ор8р)) +
— ¥ и о ~ ^орчу ~ ^оБпи^п
+В?0{ч1о ^Мп^МЫ VIЫт)
еор$р + е{Л8Бп + е^о + ^брЫ 8 {х- х{Ы,т
+
ЬорЧ1о
р
к
г
Н
Н
В
GAS INJECTION 100 MM SCF/D
OIL PRODUCTION WELL
ф H. FT. KX Ky (UNKI SW S0 .12
О 8360 FT.
0325 FT. B33S FT.
:<■ e«o FT.
1=^1 2Э45б7вЯ10
Рис. 9. Структура матрицы для сеточной области, показанной на рисунке 1, при использовании неявного метода
0*5 MJeCllON WFl I
1= У 2 J 4 5_в 7
PRODUCTION WELL
Рис. 10. Структура матрицы для сеточной области, показанной на рисунке 1, при использовании неявного метода
для нефтяной фазы;
At
Ф
Sg + RsoSo « ~г D
j, I Sg + RsoSo VXBr, + Bo
+ cgp$p + CgSwÔSw + CgSo^So f =
= V • ((Tg + Egpôp + Egs(5Sw + SSo))V$g) + V • {Tgv (dg^p)) + V • (TlgV (dgs (5SW + 5So])) +V • [(RS0 (To + EopSp + E0SwSSw + EoSoSSo) + rspTf05p)V<ï>l0] + V • [R^V (dopôp)) +
+
X I qi + v^Nw srMwj
Bg \qg + ^j=l ^m=l
RLinl + VNw VMw
B0 \q0 + 2^j=l m=l
B (Riobop + rSp) 5p
WI(jm)
egpbp + j^Sw + SSo) + egpjpj I ô(x - x
j)
j)
r (j,m)
e0$8p + e«ijSw
+ e0pbhSpbk + e,
opbh
Spbh
x x
+
С(j,m)
bgp4g
g +
Sp+
(39)
для газовой фазы.
При неявном учете скаженных источников к уравнениям (37)-(39) добавляются уравнения (33), (34) и (35). Система уравнений для блочных матриц решается с помощью библиотеки параллельных решателей LParSol версия 2.0 [12].
Расчеты тестовых задач
Далее будут представлены результаты расчетов задач из набора тестов SPE (общества инженеров нефтяников). Выполнения тестов SPE является необходимым условием при сертификации программных комплексов нефтегазовой отрасли [6].
Тест SPE1.
Постановка теста приводится в работах [6,13]. На рисунках 9 и 10 представлена моделируемая область с нанесенной на ней сеткой. На рисунке 9 приведены свойства пласта.
Коллектор представляет собой прямоугольный параллелепипед, длина которого совпадает с шириной и равна 10000 футов. Кровля расположена на глубине 8325 футов, подошва на глубине 8425 футов. Моделируемый объект состоит из трех пластов, с разными свойствами породы (см. рисунок 9). Область покрывается прямоугольной сеткой с шагом 1000 футов по осям X и Y, по оси Z размеры сеточных блоков соответствуют высотам пластов: 20, 30 и 50 футов соответственно. Более подробное описание свойств породы и флюидов представлены в [13].
Нагнетательная скважина работает в режиме заданного расхода газа с дебитом 100 • 106 фут3/сутки. Добывающая скважина работает в режиме заданного дебита равного 20000 баррелей/сутки при стандартных (поверхностных) условиях. При достижении минимального забойного давления 1000 psi, скважина переходит в режим работы с заданным забойным давлением. Период моделирования составляет 8 лет или заканчивается при достижении одного из ограничений.
На рисунках 11, 13, 15 представлены результаты сравнения пакета НИМФА с одним из КПО, представленным в этой статье. Сравнение осуществлялось по следующим параметрам: дебит нефти, давление в блоке с добывающей скважиной, давление в блоке с нагнетательной скважиной и газонефтяной фактор. На рисунках 12, 14 и 16 представлены результаты из рассматриваемой статьи, для
n
g
o
кпо
НИМФА
Рис. 11. Дебит нефти
Рис.
тьи)
6 7 8 9 10 TIME, YEARS
12. Дебит нефти (результаты из ста-
Рис. 13. Давление в ячейке с нагнетательной скважиной
■ MOBIL SHELL INTERCOMP EXXON • CMG
CMG AGREES WITH MOBIL UP TO 4 YEARS
MOBIL-INTERCOMP-
MOBIL-INTERCOMP-
10 TIME, YEARS
Рис. 14. Давление в блоке с нагнетательной скважиной (результаты из статьи)
наглядного представления количественного разброса результатов среди коммерческих программных продуктов.
Из приведенных выше графиков временных зависимостей можно сделать вывод, что пакет программ НИМФА на данном тесте получил хорошее качественное и количественное согласие с другими программными продуктами.
Тест SPE 7.
В седьмом тесте SPE рассматривается процесс разработки месторождения с помощью горизонтальных скважин [14]. Он также используется для сравнения разных подходов к моделированию горизонтальных скважин и исследованию влияния длины и дебита скважины на извлечение нефти. В своей постановке он имеет порядка 8 вариантов, в которых используются разные длины, разные дебиты и режимы работы скважин [14]. В данной статье нами будут представлены результаты только по вариантам 1а и 1Ь. Режимы работы скважин и их параметры представлены в таблице 1:
Таблица 1
Режимы работы скважин
Вариант Длины скважин (й) Дебит жидкости добывающей скважины (вТБ/ё) Давление нагнетательной скважины М
1a 900 3000 3700
1b 2100 3000 3700
Рис. 15. Давление в ячейке добывающей скважиной
MOBIL SHELL AMOCO
-INTERCOMP
- EXXON
----SSC
TIME. YEARS
Рис. 16. Давление в ячейке добывающей скважиной (результат из статьи)
Рис. 17. Газонефтяной фактор статьи)
-ЧлП —
Рис. 19. Расчетная область (вид сверху)
Рис. 20. Расчетная область (объемное изображение)
Рис. 21. Дебит нефти и суммарная нефтедобыча для варианта теста «1а»
Time (days)
Рис. 22. Дебит нефти и суммарная нефтедобыча для варианта теста «1а» (результат из статьи)
Вид расчетной области представлен на рисунках 19 и 20. Сетка имеет равномерное разбиение по оси X и сгущается по оси Y в районе скважины. Высота ячеек по оси Z соответствует высотам слоев.
Далее на рисунках 21—24 представлены результаты сравнения пакета НИМФА с одним из КПО. Сравнения осуществляется по дебиту нефти и суммарной нефтедобычи. На рисунках справа представлены результаты из статьи [13].
Анализируя полученные результаты можно сделать вывод, что на данном тесте пакет НИМФА имеет хорошее качественное и количественное согласие результатов с одним из КПО. Сравнение результатов по водонефтяному фактору, рассматриваемому в статье, не представлено, так как по этому параметру также получен хороший результат.
Основные принципы распараллеливания комплекса НИМФА
Распараллеливание комплекса НИМФА базируется на следующих принципах:
1. Используется параллелизм по данным в модели распределенной памяти;
2. Для обменов сообщениями используется интерфейс MPI;
3. Все вычислительные ядра равноправные (нет выделенного управляющего ядра);
4. Используется поячеечное наложение параобластей на один слой ячеек.
Рисунок 25 поясняет, как происходит разделение множества ячеек параобластей на подмножества.
На рисунке 25, а красной линией показана граница параобластей. На рисунке 25, б зеленым цветом показаны ячейки, лежащие на границе параобластей. Такие ячейки в комплексе НИМФА называются оболочечными. На рисунке 25, в фиолетовым цветом показаны ячейки параобласти, лежащие за границей счетных ячеек. Такие ячейки называются несчетными. Они вводятся для организации MPI-обменов. Каждая оболочечная ячейка «знает» где лежит ее образ (известен номер параобласти и
Рис. 23. Дебит нефти и суммарная нефтедобыча для варианта теста «1Ъ»
Рис. 24. Дебит нефти и суммарная нефтедобыча для варианта теста «1Ъ» (результат из статьи)
а)
б)
в)
Рис. 25. Разделение ячеек параобластей на подмножества: а) границы параобластей; б) выделение слоев ячеек, лежащих на границе; в) выделение слоев ячеек, лежащих за границами
локальный номер несчетной ячейки этой параобласти).
В комплексе НИМФА могут использоваться как явные разностные уравнения, так и неявные. В случае явных уравнений применяется метод «встречных обменов», позволяющий совмещать вычислительную работу ядра с его обменами. В случае неявных разностных схем распараллеливание происходит на уровне решения матрицы СЛАУ посредством библиотеки ЬРагёо1.
Приведем результаты эффективности распараллеливания (в режиме дробления) на тестовой задаче 8РЕ1. Количество ячеек в задаче составило порядка 1 млн. Использовалась библиотека параллельных решателей ЬРагёо1. Результаты расчетов представлены в таблице 2.
Таблица 2
Результаты расчетов задачи БРЕ1 в параллельном режиме
Количество Время счета Эффективность Ускорение
ядер (секунды) (%)
1 15005.5 100 1
96 218.7 71 69
192 122.2 64 123
384 76.3 51 197
768 40.3 48 372
На рисунке 26 приведен график зависимости эффективности от количества ядер. На рисунке 27 представлен график ускорения счета с увеличением количества ядер. Из рисунков 26,27 видно, что эффективность распараллеливания и ускорение счета в тестовой задаче 8РЕ1 являются достаточно высокими. При этом для хороших показателей эффективности необходимо обеспечить порядка 10000 ячеек на вычислительное ядро.
Рис. 26. эффективность распараллеливания тестовой задачи БРЕ1
Рис. 27. Ускорение при расчете тестовой задачи БРЕ1
Заключение
В статье были представлены основные особенности пакета программ НИМФА, которыми являются: параллельный сеточный генератор, полностью неявная схема решения уравнений модели «Black Oil» на неструктурированных сетках, высокий уровень распараллеливания вычислительных модулей. На задачах из набора тестов SPE было получено хорошее качественное и количественное согласие результатов с коммерческими программными продуктами.
Таким образом, в настоящий момент создан отечественный программный комплекс НИМФА, позволяющий проводить расчеты задач многофазной фильтрации по модели «Black Oil» в параллельном режиме.
ЛИТЕРАТУРА
1. Официальный сайт компании Schlumberger [электронный ресурс]. URL: http://www.slb.ru/sis. (дата обращения: 20.10.2014)
2. Официальный сайт компании Roxar [электронный ресурс]. URL: http://www.roxar.com. (дата обращения: 20.10.2014)
3. Официальный сайт компании Computer Modeling Group Ltd. URL: http://www.cmgroup.com/software/stars.htm (дата обращения: 20.10.2014)
4. Бутнев О. И., Пронин В. А., Сидоров М. Л., Колесников С. С., Кузнецов В. Ю. Пакет программ НИМФА-2 для решения задач многофазной фильтрации с применением суперкомпьютерных технологий // Труды XIV международной конференции «Супервычисления и математическое моделирование», РФЯЦ-ВНИИЭФ, Саров, 2012. С.112-119.
5. Сидоров М. Л., Пронин В. А. Неструктурированная призматическая дискретизация сложных геологических структур в параллельном режиме. ВАНТ. 2015.
6. Методические указания по созданию постоянно действующих геолого-технологических моделей нефтяных и газонефтяных месторождений (Часть 2. Фильтрационные модели): Утв. Заместителем Министра энергетики РФ Шелеповым В.В. 22.01.2002 г. М.: ОАО "ВНИИОЭНГ", 2003. 228 с.
7. Бутнев О. И., Бардина М. Н., Горев И. В. и др. Суперкомпьютерные технологии для нефтегазовой отрасли. Сб. статей Международной конференции RAO/CIS Offshore, г.Санкт-Петербург, 2011. С. 499-502.
8. Zgangxin Chen, Gyamrem Huan, Yuanle Ma Computational methods for multiphase flows in porous media. Dallas, Texas, Southern Methodist University, 2006.
9. Peaceman D. W. Interpretation of well-block pressure in numerical reservoir simulation // SPE Journal. 1983. V.23, № 3. P. 531-543.
10. Douglas J. Jr., Peaceman D. W., Rachford H. H. Jr., A method for calculating multi-dimensional immiscible displacement // Trans. SPE AIME 216. 1959. P. 297-306.
11. Coats K. H. Implicit compositional simulation of single-porosity and dual-porosity reservoirs, SPE 18427 // SPE Symposium on Reservoir Simulation. Houston, Texas. 1989.
12. Артемьев А. Н., Бартенев Ю. Г., Басалов В. Г., Бондаренко Ю. А., Варгин А. М., Голубев А. А., Ерзунов В. А., Ломтев А. В., Максимов А. С., Панов А. И., Прокофьев А. И., Романова М. Д., Фролова Н. В., Щаникова Е. Б. Библиотека решателей разреженных линейных систем // Труды РФЯЦ-ВНИИЭФ, Научно-исследовательское издание, Выпуск 7, Саров. 2004. C. 80-95.
13. Aziz K., Odeh A. S. Comparisn of Solution to a Three-Dimensional Black-Oil Reservoir Simulation Problem // JPT. January 1981. Vol. 33 P.13025.
14. Nghiem L. S., Collins D. A., Sharma R., Seventh SPE comparative solution project: Modeling of horizontal wells in reservoir simulation, SPE 21221 // The 11th SPE Symposium on Reservoir Simulation, Anaheim, CA. 1991.
15. Кац Р. М., Афанаскин И. В. Современные математические модели нефтяной и газовой подземной гидромеханики. Отчет. М.: НИИСИ РАН, 2012. 23 с.
