автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Методы декомпозиции для решения трехмерных задач движения жидкости в пористых средах

На правах рукописи

□□3452827

Цепаев Алексей Викторович

МЕТОДЫ ДЕКОМПОЗИЦИИ ДЛЯ РЕШЕНИЯ ТРЕХМЕРНЫХ ЗАДАЧ ДВИЖЕНИЯ ЖИДКОСТИ В ПОРИСТЫХ СРЕДАХ

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

"1

Казань-2008

003452827

Работа выполнена в Институте механики и машиностроения КазНЦ РАН

Научный руководитель: кандидат физико-математических наук,

старший научный сотрудник, Мазуров Петр Алексеевич

Официальные оппоненты:

доктор физико-математических наук, профессор Бадриев Ильдар Бурханович

доктор физико-математических наук, старший научный сотрудник, Копысов Сергей Петрович

Ведущая организация: Институт прикладной математики

им. М.В. Келдыша РАН (г.Москва)

Защита состоится «04» декабря 2008 г. в 16 часов на заседании диссертационного совета Д 212.081.21 в Казанском государственно^ университете по адресу: 420008, Казань, ул. Кремлёвская, 18, корп. 2, ауд. 218.

С диссертацией можно ознакомиться в Научной библиотеке им. Н.И. Лобачевского Казанского государственного университета

Автореферат разослан «03» ноября 2008 г.

Ученый секретарь диссертационного совета Д 212.081.21 р

д.ф.-м.н. 0<2\< Задворнов O.A.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

- разработка методов декомпозиции области для численного решения трехмерных задач фильтрации жидкости в областях сложной геометрии;

- построение алгоритмов, основанных на методах декомпозиции области, и их реализация на многопроцессорных вычислительных системах;

- численное тестирование алгоритмов с использованием различного числа процессоров при решении прикладных задач.

Научная новизна результатов.

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

- напорной однофазной фильтрации, подчиняющейся линейному закону Дарси и нелинейному закону Форхгеймера;

- напорно-безнапорной фильтрации;

- двухфазной и трехфазной фильтрации.

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

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

з

Апробация работы. Основные результаты диссертационной работы докладывались и обсуждались на XI международной конференции по вычислительной механике и современным прикладным программным системам (ВМСППС'2001, Москва-Истра, 2001), IV научно-практической конференции молодых ученых и специалистов Республики Татарстан (Казань, 2001), VIII Четаевской международной конференции «Аналитическая механика, устойчивость и управление движением» (Казань,

2002), конференции «Современные проблемы гидрогеологии и гидромеханики» (Санкт-Петербург, 2002), Всероссийской конференции «Высокопроизводительные вычисления и технологии (ВВТ-2003)» (Ижевск,

2003), XVII сессии международной школы по моделям механики сплошной среды (Казань, 2004), международном семинаре «Супервычисления и математическое моделирование» (Саров, 2006), итоговых научных конференциях КазНЦ РАН, на семинарах Института механики и машиностроения КазНЦ РАН.

Работа выполнена в рамках программы фундаментальных исследований Президиума РАН "Параллельные вычисления на многопроцессорных вычислительных системах".

Структура диссертационной работы. Диссертация состоит из введения, пяти разделов, выводов и списка литературы. Общий объем диссертации составляет 97 страниц, включая 18 таблиц и 23 рисунков.

Публикации. По теме диссертации опубликовано 12 работ, список которых приведен в конце автореферата.

СТРУКТУРА И КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ.

Во введении дается обзор литературы, формулируются цели исследования и положения, выносимые на защиту.

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

В п. 1.1 даётся постановка задачи. Необходимо определить поле напора h в ограниченной области D из решения уравнения

divKgradh = 0 (1.1)

при граничных условиях

h = Н0 на Г,, (1.2)

-Kah/3n = qr„ на Гг, (1.3)

- jKdh/dnds = Qk при h|Si = const или h|Si=Hk, k=l,...,N, (1.4)

где Г = Г, иГ2 - внешняя граница области D, Г, пГ; =0, К - коэффициент фильтрации, Sk - суммарная поверхность интервалов вскрытия к-ой скважины, Но и Нк - заданные напоры, q,.n - заданное значение нормальной

составляющей скорости фильтрации, 0к - расход к-ой скважины, N -число скважин. Уравнения (1.1)-(1.4) описывают напорную однофазную стационарную фильтрацию, подчиняющуюся закону Дарси.

В п. 1.2 строится разбиение области решения задачи на подобласти.

_ _ N _

Для области Б имеем следующие соотношения: и(иОк),

\ , Б, п = ук, где Ок - прискважинные подобласти. Определим области

N

Ук как области интервалов вскрытия скважин. Объединение иУк является

дополнением многосвязной области Б до односвязной. Поверхности интервалов вскрытия являются внутренними границами области О. Весь пласт, как односвязная область, покрывается грубой сеткой. От узлов грубой сетки, расположенных на граничной поверхности прискважинных подобластей Ук> строятся сетки, сгущающиеся к интервалам вскрытия скважин. Грубая сетка вне сгущающихся участков является основной грубой сеткой, а грубые сетки со, на сгущающихся участках считаются дополнительными грубыми сетками. Грубая сетка вне сгущающихся участков и сетки, сгущающиеся к интервалам вскрытия скважин, образуют сетку О, в узлах которой требуется определить сеточную функцию Ьа из решения системы (1.1)-(1.4). Без декомпозиции области решение Ь определяется на сетке £~2. В предлагаемом методе с декомпозицией области решение Ь в области представляется как Ь,, а в областях Ок представляется в виде суммы двух решений: одно Ь,к определено в узлах сгущающегося участка сетки, другое Ь31с - в узлах дополнительной грубой сетки сок в области Эк и Ук.

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

В п. 1.4 строится алгоритм решения задачи в общем случае на основе предлагаемого метода декомпозиции области. Даются постановки задач для И,, Ь2к, Ь.к. Решения Ь,, Ь2к и Ьзк при известных граничных значениях Ь, , Ьг , на независимо определяются из соответствующих систем уравнений. При выполнении условий на границах раздела ук относительно напоров ЬГ1 = ЬГг1 +ЬГ и относительно нормальных составляющих скоростей фильтрации я1кл + ч2к11 + я1Ь1 = 0 для определения решений Ьр Ь2к, Ь,к достаточно задания граничных значений . При значениях Ь. = О

система уравнений эквивалентна исходной системе уравнений. Для решения соответствующих систем уравнений строится итерационный процесс. С увеличением номера итерации 1 значения Ь|,к стремятся к нулю, и ¡-ое приближение напора, определяемое решением Ь] во внескважинной области

и решениями в прискважинных подобластях, стремится к решению системы (1Л)-(1.4).

Данный алгоритм для определения напора можно записать в виде:

1) В начальном приближении Ь°к=0, а Ь° и совместно определяются на грубой и дополнительных грубых сетках из решения системы (1.1)-(1.4).

2) На ¡-ой итерации независимо решаются системы уравнений на сгущающихся участках для определения Ь'2к при граничных условиях

на границах ук с учетом условия (1.4).

3) Совместно определяются решения Ь|, Ь'зк на грубой и дополнительных грубых сетках с условиями я'1кп + q'21cn + = 0 и

= Ь^ + Ь'г<. на границах раздела ук без учета (1.4).

В п. 1.5 приводятся результаты численных экспериментов. Рассматривался пятислойный пласт с различными толщинами слоев и коэффициентами фильтрации. Пласт считался . непроницаемым за исключением кровли, на которой задавалось граничное условие 2-го рода qI,=3.4x 10"4м/сут, и противоположных участков боковой поверхности пятого слоя с граничными условиями 1-го рода Н0 = 10 м и Н„ = 40 м. Для аппроксимации уравнений использовался метод Галеркина, полученная система алгебраических уравнений решалась методом сопряженных градиентов с предобуславливающей матрицей, построенной с помощью неполного разложения Холесского. Задачи решались на многопроцессорной вычислительной системе МВС-1000. В табл. 1 приведено время решения задач с различным числом скважин, полученное алгоритмом с декомпозицией области (различное число процессоров) и без декомпозиции области (один процессор). В обоих случаях под решением понимается определение сеточной функции Ь0 с одной и той же заданной точностью.

Таблица 1. Время решения задач при заданных расходах на скважинах.

Число Число Число Решение задачи с Решение задачи

скважин узлов процессоров декомпозицией области без декомпозиции области

1 15978 1 17сек 16сек

10 46866 1 1мин 21 сек 1мин 25сек

5 55сек

50 184146 1 4мин 57сек 9мин 26сек

5 2мин 05сек

14 1мин 5 Осек

100 355746 1 12мин 45 сек 176мин Збсек

5 5мин 14сек

14 2мин 25сек

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

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

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

двучленным нелинейным законом Форхгеймера £гас! Ь = -(А + Ву)у, где А=1/К, К - коэффициент фильтрации, V - модуль вектора скорости фильтрации V, В - константа пористой среды. Задача решается в следующей постановке: требуется определить поле напора Ь в области О из решения системы уравнений

(Пуу = ОВО, (2.1)

ЕгаёЬ = -(А + Ву)у в в,, к = 1,...,М, (2.2)

К^ас1 Ь = -у в Э \ ^ (2.3)

при граничных условиях

Ь = Н0 на Г,, (2.4)

У„=УГ»НаГ2' (2-5)

]упё5 = (Зк при И 1^ = 001^ или Ь|8к = Нк, к = 1,...,Ы, (2.6)

«I

где йк ={хе 11е(у)> Яе^} - области выполнения нелинейного закона фильтрации, Г = Г,иГ2 - внешняя граница области В, Г, п Г, =0, -суммарная поверхность интервалов вскрытия к-ой скважины, Н0 и Нк -заданные напоры, V - скорость фильтрации, уо - нормальная составляющая скорости фильтрации, уП1 - заданное значение нормальной составляющей скорости фильтрации, <Зк - расход к-ой скважины, N - число скважин. Границы областей Ск заранее неизвестны и должны определяться в процессе решения. В постановке задачи предполагается йк с Эк. В случае нарушения этого условия границы прискважинных зон Эк могут быть расширены вплоть до общих границ (В„=0).

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

применения нелинейного закона фильтрации берётся с предыдущего шага. Затем формируются области которые состоят из конечных элементов, где числа Рейнольдса Яе(у) больше критического значения Ле^. Итерационный процесс останавливается в случае выполнения условия в!"' = и достижения заданной точности решения системы уравнений.

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

В п. 2.4 приводятся результаты численных экспериментов. Рассматривался пятислойный пласт с различными толщинами слоев и коэффициентами фильтрации. Задача решалась с использованием закона Форхгеймера при различных критических числах Рейнольдса с заданными на скважинах расходами и напорами. В табл. 2 приведено время решения задач фильтрации при заданных расходах на скважинах. Расчеты проводились на многопроцессорной вычислительной системе МВС-1000.

Таблица 2. Время решения нелинейных задач фильтрации жидкости при заданных расходах на скважинах при 11екр = 0.01.

Число скважин Число узлов Число процессоров Решение задачи с декомпозицией области Решение задачи без декомпозиции области

1 15978 1 40сек 2мин

1 5мин 18мин

10 46866 5 1мин ЗОсек

10 1мин

1 28мин ЮОмин

50 184146 5 бмин

10 4мин

1 85мин 440мин

100 355746 5 21 мин

10 13 мин

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

В п.3.1 дается постановка задачи. Задача решается в области ОсБ. Область О представляет собой пласт, ограниченный кровлей, подошвой, боковыми поверхностями и поверхностями интервалов вскрытия скважин У^. На участках Г, и Гг границы области Ц задаются граничные условия первого и второго рода. Область решения О определяется в процессе решения с учетом депрессионных поверхностей в прискважинных зонах. Участки граничных условий первого и второго рода Г0], Г02 определяются положением верхней границы Гн6 области в, при этом Гв5 сГы, Гн5 = Гя + Г6, где Ь > ъ на Г„ и Ь = ъ на Г5. Решается следующая задача: область б и поле напора Ь в в определяются из решения уравнения

сПу1^гасШ = 0 в в (3.1)

при граничных условиях

Ь = ЬГ на Г01, (3.2)

-К5Ь/еп = Чг„ наГС2, (3.3)

Ь>2 наГ„,Ь = 2наГ6,Г„+Г6=Гк5сГ02, (3.4)

- |КдЬ/дпсЬ = С>к при Ь|51 = соп51 или ИНк, к=1,...,Ы, (3.5)

где К - коэффициент фильтрации, Ьг и Нк - заданные напоры, - заданное значение нормальной составляющей скорости фильтрации, Бк - суммарная поверхность интервалов вскрытия к-ой скважины, СЬ - расход к-ой скважины, N - число скважин.

В п. 3.2 строится итерационный алгоритм для решения задачи напорно-безнапорной фильтрации жидкости, основанный на методе декомпозиции области, аналогичном изложенному в п. 1.4. В каждой прискважинной подобласти решение представляется в виде суммы двух решений. Одно решение определено на сгущающейся сетке, другое - на дополнительных грубых сетках. Независимо решаются системы уравнений для сгущающихся участков сетки, затем совместно решается система уравнений для основной и дополнительных грубых сеток. При решении уравнений для сгущающихся участков сетки определяются части верхней поверхности, на которой Ь < г, уточняются координаты ъ депрессионной поверхности из условия ъ = Ь, и все сетки перестраиваются от депрессионной поверхности. На дополнительных грубых сетках режим фильтрации считается напорным.

В п. 3.3 приводятся результаты численных экспериментов. Задачи решались на многопроцессорной вычислительной системе МКВС-Е112. Результаты расчетов показали, что с возрастанием числа сгущающихся участков сетки увеличивается разница во времени решения задачи по алгоритму без декомпозиции области и по предложенному алгоритму с декомпозицией области даже на однопроцессорном компьютере. На рис. 1 приведены графики, показывающие зависимость времени решения задачи от числа процессоров и числа сгущающихся участков сетки.

▲ - решение задачи по предложенному алгоритму с числом процессоров, равным числу сгущающихся участков • - решение задачи по предложенному алгоритму на одном процессоре ■ - решение задачи без декомпозиции области пуч - число сгущающихся участков сетки 1 - время решения

Рисунок 1. Расчетное время решения задачи.

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

В п. 4.1 описывается схема метода решения уравнения насыщенности, основанного на декомпозиции области. Процесс построения алгоритма для определения насыщенности на п-ом временном шаге состоит из следующих этапов:

1) По явной схеме вычисляются фазовые расходы, выходящие из ячеек грубой сетки, окружающих ячейки сгущающихся участков прискважинных зон.

2) Определяются насыщенности по неявной схеме для ячеек сгущающихся участков, при этом фазовые расходы, входящие в прискважинные зоны, являются граничными условиями при решении указанных систем.

3) По неявной схеме вычисляются фазовые расходы, выходящие из ячеек сгущающихся участков в грубые ячейки.

4) Подправляются насыщенности для ячеек грубой сетки.

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

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

сНу((К„ + К^гас1р) = 0, (4.1)

сИу(ч,,)+т08„/а = О (4.2)

при граничных условиях

Р = Рг на Г,, (4.3)

-(ко + к,„)ар/ап = яГп на г2, (4.4)

ри.= Рк, к = 1,...,н, (4.5)

= Б, на Г, № \*г 3 (4-6)

и начальном условии

= в Б, (4.7)

где р = р(х,у,г) - давление, \'г - вектор скорости фильтрации вытесняющей жидкости (воды), q =-(Ко + Kl>)gradp - вектор скорости фильтрации, -нефтенасыщенность, Б,,, - водонасыщенность, 8о + = 1, Ко =Ко(8„) = кГо/цо, = = к^ - фазовые подвижности,

^ =^(5о)- относительные фазовые проницаемости, к -абсолютная проницаемость, - динамические вязкости фаз, т —

пористость, Г,иГ2 = Г - внешняя граничная поверхность области Б, Г, -часть поверхности Г, через которую жидкость поступает в пласт, 5Ук -поверхность интервала вскрытия к-ой скважины, Рк - заданное давление на к-ой скважине, - заданное значение нормальной составляющей вектора скорости фильтрации, N - число добывающих скважин. На п-ом временном шаге поле давлений рг вычислялось по методу декомпозиции области,

и

описанному в п. 1.4. Для нахождения насыщенности применялся метод декомпозиции области, описанный в п. 4.1.

В п. 4.3 приводятся результаты численных экспериментов. Брался десятислойный пласт с различными толщинами слоев и абсолютными проницаемостями. Кровля пласта считалась непроницаемой, на боковых поверхностях, подошве пласта и на скважинах задавалось давление. Начальная насыщенность Б,, = 1, на боковой поверхности 8о = 1, на подошве = 1. Относительные фазовые проницаемости брались линейными функциями от насыщенностей.

Аппроксимация систем уравнений проводилась методом контрольных объемов. При решении системы алгебраических уравнений для определения поля давлений использовался метод сопряженных градиентов с предобуславливающей матрицей, построенной с ¡помощью неполного разложения Холесского. При решении системы алгебраических уравнений для определения насыщенности для ячеек сгущающихся участков сетки использовался метод Зейделя. !

Задачи решались на многопроцессорной вычислительной системе МКВС-Е112. На рис. 2 приведены результаты решения задачи без декомпозиции области на одном процессоре и с декомпозицией области по предложенному алгоритму. Показана эффективность алгоритма при решении задач с большим числом сгущающихся участков сетки по сравнению с

А - решение задачи по предложенному алгоритму с числом процессоров, равным числу сгущающихся участков • - решение задачи по предложенному алгоритму на одном процессоре ■ - решение задачи без декомпозиции области пуч - число сгущающихся участков сетки I — время решения

Рисунок 2. Расчетное время решения с декомпозицией области по давлению и насыщенности.

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

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

решением без декомпозиции области.

1 50 100

закону Дарси. Считается, что нефть и вода несжимаемы, отсутствует массообмен между нефтяной и газовой фазами. Гравитационные и капиллярные силы не учитываются. Насыщенности в пласте в начальный момент времени считаются известными. На внешней поверхности пласта задаются граничные условия 1-го или 2-го рода для давления и известные насыщенности фаз на участках внешней поверхности, через которые поступают флюиды. На скважинах задаётся забойное давление.

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

В п. 5.3 приводятся результаты численных экспериментов. Задачи тестировались на многопроцессорной вычислительной системе МКВС-Е112. На рис. 3 приведены результаты решения в случае, когда задача решалась с декомпозицией области по предложенному алгоритму и без декомпозиции области. Показана эффективность предложенного алгоритма при решении задач с большим числом сгущающихся участков сетки по сравнению с решением без декомпозиции.

А - решение задачи по предложенному алгоритму с числом процессоров, равным числу сгущающихся участков • - решение задачи по предложенному алгоритму на одном процессоре ■ - решение задачи без декомпозиции области пу,, - число сгущающихся участков сетки 1: — время решения

Рисунок 3. Расчетное время решения с декомпозицией области по давлению и насыщенностям.

В заключении приводятся основные результаты работы.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ

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

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

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

СПИСОК ОПУБЛИКОВАННЫХ РАБОТ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Цепаев A.B. Один из алгоритмов метода разделения области для решения уравнения фильтрации с использованием принципа суперпозиции // Тезисы докладов IV научно-практической конференции молодых ученых и специалистов Республики Татарстан. - Казань, 2001. - С. 55

2. Цепаев A.B. К решению задач фильтрации несжимаемой жидкости в трехмерных пластах с гидродинамически несовершенными скважинами / П.А. Мазуров, A.B. Цепаев // Тезисы международной конференции по вычислительной механики и современным прикладным программным системам. - Москва-Истра, 2001. - С. 251.

3. Цепаев A.B. К решению задач фильтрации несжимаемой жидкости в трехмерных пластах с гидродинамически несовершенными скважинами / П.А. Мазуров, A.B. Цепаев // Математическое моделирование. - 2002. - Т. 14. -№9. -С. 121-123.

4. Цепаев A.B. Метод суперпозиции для решения задач фильтрации жидкости в трехмерных пластах с гидродинамически несовершенными скважинами / П.А. Мазуров, A.B. Цепаев // Современные проблемы гидрогеологии и гидрогеомеханики. Сб. докл. конф. - СПб., 2002. - С. 471476.

5. Цепаев A.B. Использование принципа суперпозиции и радиальных потоков в решении трехмерных задач фильтрации с большим числом скважин / П.А. Мазуров, A.B. Цепаев // Тезисы международной Четаевской конференции «Аналитическая механика, устойчивость и управление движением». - Казань, 2002. - С.269.

6. Цепаев A.B. Решение трехмерных задач фильтрации жидкости с большим числом гидродинамически несовершенных скважин на многопроцессорных вычислительных системах / П.А. Мазуров, A.B. Цепаев // Тезисы

Всероссийской конференции «Высокопроизводительные вычисления и технологии». - Ижевск, 2003. - С. 55

7. Цепаев A.B. Метод решения нелинейных задач фильтрации жидкости в трехмерных пластах с гидродинамически несовершенными скважинами / П.А. Мазуров, A.B. Цепаев // Математическое моделирование. - 2004. - Т.16. -№3. - С. 33-42.

8. Цепаев A.B. Решение трехмерных задач фильтрации жидкости на МВС-1000/16 на сетках со сгущающимися участками / П.А. Мазуров, A.B. Цепаев // Актуальные проблемы механики сплошной среды. - Казань, 2004. - С 45-56.

9. Цепаев A.B. Алгоритм решения трехмерных задач напорно-безнапорной стационарной фильтрации жидкости со сгущающимися участками сетки / Д.А. Губайдуллин, П.А. Мазуров, A.B. Цепаев // Вычислительные методы и программирование. - 2005. - Т.6. - №2. - С. 217 - 225.

10. Цепаев A.B. Алгоритмы для распараллеливания решения задач двухфазной фильтрации жидкости на сетках со сгущающимися участками / П.А. Мазуров, A.B. Цепаев // Вычислительные методы и программирование. - 2006. - Т.7. - №2. - С. 115 - 123.

11. Цепаев A.B. Алгоритмы для распараллеливания решения задач двухфазной фильтрации жидкости на сетках со сгущающимися участками / П.А. Мазуров, A.B. Цепаев // Тезисы международного семинара «Супервычисления и математическое моделирование». - Саров, 2006. -С. 56.

12. Цепаев A.B. Алгоритмы распараллеливания на сгущающихся сетках в задачах трехфазной фильтрации жидкости / Д.А. Губайдуллин, А.И. Никифоров, A.B. Цепаев // Вычислительные методы и программирование. -2007. - Т.8. - №2. - С. 360 - 366.

Отпечатано с готового оригинал-макета в типографии Издательства Казанского государственного университета Тираж 130 экз. Заказ 103/10

420008, ул. Профессора Нужина, 1/37 тел.: 231-53-59, 292-65-60

Основные обозначения.

Введение.

Раздел 1. Задача однофазной напорной фильтрации жидкости в трехмерном пласте, подчиняющейся закону Дарси.

1.1. Постановка задачи.

1.2. Разбиение области решения на подобласти.

1.3. Описание метода декомпозиции области для решения уравнения по напору.

1.4. Алгоритм решения задачи с декомпозицией области.

1.5. Численные эксперименты.

1.6. Выводы. 39 Раздел 2. Задача однофазной напорной фильтрации жидкости, подчиняющейся закону Форхгеймера.

2.1. Постановка задачи.

2.2. Алгоритм решения задачи без декомпозиции области.

2.3. Алгоритм решения задачи с декомпозицией области.

2.4. Численные эксперименты.

2.5. Выводы. 52 Раздел 3. Задача однофазной напорно-безнапорной фильтрации жидкости в трехмерном пласте.

3.1. Постановка задачи.

3.2. Алгоритм решения задачи с декомпозицией области.

3.3. Численные эксперименты.

3.4. Выводы. 63 Раздел 4. Задача двухфазной фильтрации жидкости в трехмерном пласте.

4.1. Описание метода декомпозиции области для решения уравнения по насыщенности.

4.2. Постановка и алгоритм решения задачи двухфазной фильтрации.

4.3. Численные эксперименты.

4.4. Выводы. 78 Раздел 5. Задача трехфазной напорной фильтрации жидкости в трехмерном пласте.

5.1. Постановка задачи.

5.2. Алгоритм решения задачи с декомпозицией области.

5.3. Численные эксперименты.

5.4. Выводы. 88 Заключение. 89 Литература.

Основные обозначения.

D - область решения

Dk - прискважинная область k-ой скважины D0=D\(uDk) - внескважинная область

GkcDk - область выполнения нелинейного закона фильтрации в к-ой прискважинной зоне G - область решения (раздел 3) yk = Dk n D0- часть граничной поверхности к-ой прискважинной зоны Г - граница области D

Г, - участок граничной поверхности Г, на которой задаются граничные условия 1-го рода

Г2 - участок граничной поверхности Г, на которой задаются граничные условия 2-го рода

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

Г02 - участок граничной поверхности области решения G, на котором задаются граничные условия второго рода (раздел 3)

Гн6 - верхняя граничная поверхность области G

Гн - часть поверхности Гиб, соответствующая напорному участку

Гб - часть поверхности Гн6, соответствующая безнапорному участку

Vk - область интервалов вскрытия k-ой скважины

Sk - поверхность интервала вскрытия (раздел 1, 2, 3)

5Vk - поверхность интервала вскрытия (раздел 4, 5)

Q - основная сетка со сгущающимися участками go k - дополнительная грубая сетка k-ой прискважинной зоны h - напор h,, h2k, h3k - решения на основной грубой сетке во внескважинной области, на k-ом сгущающемся участке и на k-ой дополнительной грубой сетке соответственно

Нк - заданный напор (раздел 1, 2, 3) hr t, ЬГ к - значения решений h,, h2k, h3k на поверхности ук р - давление

Qk - расход к-ой скважины п - внешняя нормаль к поверхности области D q - вектор скорости фильтрации (разделы 1, 3, 4, 5) qn - нормальная составляющая скорости фильтрации qrn - заданное значение нормальной составляющей скорости фильтрации v - вектор скорости фильтрации (раздел 2) vn - нормальная составляющая скорости фильтрации vrn - заданное значение нормальной составляющей скорости фильтрации

К - коэффициент фильтрации

Kw, KD, Kg - фазовые подвижности воды, нефти и газа соответственно

Sw,S0,Sg- насыщенности воды, нефти и газа соответственно fw, fo, f - относительные фазовые проницаемости воды, нефти и газа соответственно

M'w' Mo' M'g " динамические вязкости воды, нефти и газа соответственно m - пористость

Re (v) - число Рейнольдса

С развитием компьютерной техники [8] появилась возможность использования многопроцессорных вычислительных систем для решения сложных математических задач. К таким задачам относятся, в частности, задачи фильтрации жидкости в нефтяных и водоносных пластах вскрытых системой скважин. При численном решении этих задач одной из основных трудностей является решение систем алгебраических уравнений большой размерности, что обусловлено трехмерностью объекта и необходимостью сгущения сетки в прискважинных зонах. Размерности сгущающихся сеток при этом часто достигают размерности грубой сетки. Одним из быстро развивающихся методов решения таких задач являются методы декомпозиции, то есть разделение задачи на части, которые могут быть параллельно решены на нескольких процессорах. С помощью декомпозиции можно сократить общее время счета задачи, и более эффективно использовать доступную оперативную память. Основным подходом в создании таких алгоритмов являются методы декомпозиции области: метод подструктур, методы Шварца, поэлементная декомпозиция, много фронтальная декомпозиция и т.д. [32, 46, 50, 53-56, 58, 62, 73, 74, 76].

В основе метода подструктур [73] лежит идея независимого расчета отдельных подобластей и последующего учета их взаимодействия. В этом методе исходная область разбивается на подобласти, имеющие лишь общую границу. Метод подструктур является прямым методом. В каждой из подобластей общее решение записывается в виде суперпозиции двух решений. Первое определяется из решения системы уравнений, соответствующей каждой из подобластей с нулевыми граничными условиями с учетом заданных источников внутри подобласти, второе определяется исходя из суперпозиции и зависит от вектора граничных значений, полученного из системы уравнений с так называемой матрицей граничных значений (или дополнение Шура) [48]. Матрица граничных значений имеет меньшее число обусловленности, чем матрица исходной системы уравнений, однако, число обусловленности может быть достаточно большим, чтобы решить систему за приемлемое время.

Одними из наиболее популярных методов разделения области с перекрывающимися областями являются методы Шварца (итерационные методы подструктур). Альтернирующий метод Шварца известен в вычислительном математике уже многие годы [21, 34, 36, 37]. Он всегда рассматривался, прежде всего, как метод, позволяющий сводить решение исходной задачи в области со сложной формой границы к последовательности задач в подобластях, форма которых достаточно простая. В настоящее время разработаны достаточно эффективные алгоритмы численного решения ряда задач математической физики именно для случая простых областей. С учетом этого обстоятельства в последние годы осуществляется поиск модификаций методов Шварца, обладающих более высокой скоростью сходимости по сравнению с их классическим вариантом [7, 13, 16, 19, 21, 57]. Если рассматривать область решения D как объединение двух перекрывающихся подобластей D,, D2, и решать систему дифференциальных уравнений в частных производных

Lu = f в D, u = g на Г = 5D, то каждая итерация альтернирующего метода Шварца состоит в последовательном решении подзадач для и,

Lu" =f в D,, и" = g на dD, \ Г,, и; = и"4 на Г,,

И для и2

Lu2 = f в D2, и2=§наШ2\Г2, и2 = и" на Г2. где Г; - часть границы dDi принадлежащей области Dj? i Ф j.

Для многих задач возникает ситуация, когда строятся сетки отдельно для каждой подобласти, и в зоне пересечения сетки не совпадают. Для таких случаев используется мультипликативный метод Шварца [55]. При решении системы линейных уравнений Au = f мультипликативный метод Шварца записывается в виде un+,u =un +

О О f - Au"),

Un+I = un+U2 +

О 0 ^ 0 f - Au ).

В этих методах каждый итерационный шаг состоит из последовательного вычисления функции и в подобластях, что не позволяет распараллеливать задачу при решении на многопроцессорных машинах. Аддитивный метод Шварца [34] может рассматриваться как параллельная версия мультипликативного метода Шварца и записывается в виде un+1 = и" + (В, + В2 )(f - Au"), где В, =RjrA~1Ri, R; - оператор проекции во внутренние узлы подобласти D;.

В данном методе системы уравнений для подобластей решаются независимо, поэтому могут решаться параллельно. Аддитивный метод Шварца может рассматриваться как обобщение блочного метода Якоби, а мультипликативный - блочного метода Гаусса-Зейделя.

Другим методом декомпозиции области является поэлементная декомпозиция [18]. Одним из главных достоинств итерационных методов в сочетании с методом конечных элементов является то, что можно не строить глобальную матрицу системы уравнений в явном виде. Глобальная матрица может быть записана в блочно-диагональном виде. Такую матрицу можно хранить отдельно на каждом процессоре, вводя матрицу связности, которая отражает зависимость между элементами и узлами для всей области в целом. В работе [72] показано как реализуется скалярное произведение, а также вычисление произведения глобальной матрицы на вектор, участвующих в итерационных методах, без построения глобальной матрицы. Данный метод является методом поэлементной декомпозиции и эффективно реализуется на параллельных машинах с распределенной памятью [72].

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

1) меньшие затраты времени решения и памяти в случае, когда матрица коэффициентов имеет ленточный характер,

2) произвольность способа нумерации узлов,

3) возможность производить локально сгущение сетки без перенумерации,

4) в отличие от других прямых методов многофронтальный метод перспективен для распараллеливания.

В первых исследованиях по методам разделения области на подобласти, имеющих лишь общую границу, была отмечена необходимость введения в итерационный процесс, лежащий в основе метода, некоторых параметров. Эти параметры необходимо было выбирать так, чтобы процесс сходился. Различные подходы по выбору оптимальных значений этих параметров предложены В.Э. Канцельсоном, В.В. Меньшиковым [16], Л.Б. Цвиком [40]. Параметры итерационного процесса выбирались на основе одношаговой минимизации квадратичного функционала. В работах М.Е. Дмитриенко, Л.А. Оганесяна [13], В.Г. Осмоловского, В.Я. Ривкинда, A.M. Мацокина [33] рассматривались стационарные итерационные методы разделения области в применении к дифференциальным уравнениям и их разностным аналогам. В работах В.И. Лебедева, В.И. Агошкова [1] исследованы эффективные нестационарные итерационные алгоритмы разделения области с переменными параметрами, найдены оптимальные наборы этих параметров. В.И. Лебедевым, В.И. Агошковым введен специальный класс операторов - операторов Пуанкаре -Стеклова, изучен ряд свойств этих операторов. На их основе разработана общая методология конструирования и исследования метода разделения области. В.В. Смеловым и В.И. Агошковым методы разделения области изучались в применении к задачам теории переноса излучения. Исследованием данных методов в нестационарных задачах посвящены работы Ю.А. Кузнецова [19], Г.И. Марчука, С.Н. Булеева, В.И. Агошкова и др. [1, 21, 31]

К прямым методам разделения области на подобласти, имеющим лишь общую границу, также можно отнести balancing domain decomposition method (BDD) и balancing domain decomposition by constraints (BDDC) [47, 65-70]. BDD, BDCC - это итерационные методы для нахождения решения системы линейных алгебраических уравнений с симметричной, положительно определенной матрицей, построенной методом конечных элементов. В двойственных методах, таких как finite element tearing and interconnect (FETI) [43, 45, 51, 52], непрерывность решений вдоль границ соседних подобластей достигается за счет введения лагранжевых множителей [43]. Метод FETI-DP является сочетанием прямых и двойственных методов. Данный метод был предложен Charbel Farhat [51, 52], а впервые математический анализ был проведен Mandel J. и Tezaur R.

В данной работе предлагаются новые алгоритмы решения различных задач фильтрации жидкости в пласте с гидродинамически несовершенными скважинами, основанные на методах декомпозиции области. В каждой из задач выделяются прискважинные подобласти, имеющие между собой нулевое пересечение. Весь пласт, как односвязная область, покрывается грубой сеткой. От узлов грубой сетки, расположенных на граничной поверхности прискважинных подобластей, строятся сетки, сгущающиеся к интервалам вскрытия скважин. Грубая сетка вне сгущающихся участков и сетки, сгущающиеся к интервалам вскрытия скважин, образуют основную сетку. Для решения линейной задачи напорной фильтрации жидкости предлагается метод декомпозиции области, основанный на согласовании независимых решений систем уравнений на сгущающихся участках сетки за счёт введения дополнительных грубых сеток. Предложенный метод также используется при решении задач с нелинейным законом фильтрации, в задачах напорно-безнапорной фильтрации, двухфазных и трехфазных течений. Постановки задач взяты из [3, 4, 17, 31, 35]. Для определения поля насыщенности предлагается метод декомпозиции области, основанный на согласовании независимых решений уравнений на сгущающихся участках сетки за счет введения элементов явной и неявной схем.

Задачи, описанные в разделах 1 - 3, аппроксимируются методом конечных элементов [9, 14, 38, 39], а в разделах 4, 5 - методом контрольных объемов [59, 66, 71, 75]. Системы линейных алгебраических уравнений, полученные в результате аппроксимации исходных дифференциальных уравнений, можно решать прямыми или итерационными методами. На практике при большом числе переменных обычно используются итерационные методы. Одним из эффективных методов решения систем является метод сопряженных градиентов с различными предобуславливателями [60]. Скорость сходимости метода сопряженных градиентов можно увеличить за счет выбора предобуславливающей матрицы. Однако, построение эффективной предобуславливающей матрицы ведет, как правило, к большему расходу оперативной памяти [60]. В [60] при решении линейных и нелинейных задач фильтрации проводится сравнение метода сопряженных градиентов с тремя различными предобуславливающими матрицами, построенными методом неполного разложения Холесского, модифицированным методом неполного разложения Холесского и полиномиальным методом. В [64] для построения предобуславливающей матрицы используются диагональное масштабирование, неполное разложение Холесского, неполная факторизация, модифицированная неполная факторизация. В данной работе для решения систем линейных алгебраических уравнений с симметричной матрицей использовался метод сопряженных градиентов с предобуславливающей матрицей, построенной неполным разложением Холесского. В задачах по определению поля насыщенности матрица системы уравнений несимметричная, для решения таких систем использовался метод Зейделя [5].

Цели диссертационной работы:

- разработка методов декомпозиции области для численного решения трехмерных задач фильтрации жидкости в областях сложной геометрии;

- построение алгоритмов, основанных на методах декомпозиции области, и их реализация на многопроцессорных вычислительных системах;

- численное тестирование алгоритмов с использованием различного числа процессоров при решении прикладных задач.

Структура и краткое содержание работы.

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

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

В п. 1.2 строится разбиение области решения задачи на подобласти. Весь пласт, как односвязная область, покрывается грубой сеткой. От узлов грубой сетки, расположенных на граничной поверхности прискважинных подобластей, строятся сетки, сгущающиеся к интервалам вскрытия скважин. Грубая сетка вне сгущающихся участков и сетки, сгущающиеся к интервалам вскрытия скважин, образуют общую сетку Q, в узлах которой требуется определить сеточную функцию hn. Грубая сетка вне сгущающихся участков является основной грубой сеткой, а грубые сетки сок на сгущающихся участках считаются дополнительными грубыми сетками. Область решения рассматривается как многосвязная область, поверхности интервалов вскрытия скважин являются внутренними границами области. Для каждой прискважинной зоны решение h представляется в виде суммы двух решений: одно h2k определено в узлах сгущающегося участка сетки, другое h3k - в узлах дополнительной грубой сетки, к - номер сгущающегося участка.

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

В п. 1.4 строится алгоритм решения задачи в общем случае на основе предлагаемого метода декомпозиции области. Даются постановки задач для h,, h2k, h3k, где h, - решение на основной грубой сетке, h2k - решения на сгущающихся участках, h3k - решения на дополнительных грубых сетках. Решения h,, h2k и h3k при известных граничных значениях hr , hr , hr на границах раздела ук прискважинных подобластей и внескважинной области независимо определяются из соответствующих систем уравнений. При выполнении условий на границах раздела относительно напоров hri = hr t + hr и относительно нормальных составляющих скоростей фильтрации Япш +Я2кп +Чзкп = О Д3151 определения решений h,, h2k, h3k достаточно задания граничных значений hrji. При значениях = 0 система уравнений эквивалентна исходной системе уравнений. Для решения этих систем уравнений строится итерационный процесс. На каждой итерации граничные значения h^ сносятся на граничные значения h^. В результате h^ с ростом i стремятся к нулю, и i-oe приближение напора определяется решением h) во внескважинной области и решениями h'2k в прискважинных подобластях.

Рассматривается критерий прерывания решений задач для h2k и совместного решения задач для hj, h!,k, основанный на поведении дисбаланса расходов жидкости в узлах границ раздела ук.

В п. 1.5 приводятся результаты численных экспериментов. Для аппроксимации алгебраических уравнений использовался метод Галеркина, для решения сеточной системы алгебраических уравнений использовался метод сопряженных градиентов с предобуславливателем, для построения которого использовался метод неполного разложения Холесского. Рассматривался пятислойный пласт с различными толщинами слоев и коэффициентами фильтрации. Пласт считался непроницаемым за исключением кровли, на которой задавалось граничное условие 2-го рода, и противоположных участков боковой поверхности пятого слоя с граничными условиями 1-го рода. Рассматривалось два типа решений: решение системы без декомпозиции области (один процессор) и решение с декомпозицией области (различное число процессоров). Показана эффективность алгоритма при решении задач с большим числом сгущающихся участков сетки по сравнению с решением без декомпозиции области. Алгоритм апробирован на многопроцессорной вычислительной системе МВС-1000 с различным числом процессоров и гидродинамически несовершенных скважин.

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

В п. 2.1 дается постановка задачи. Пласт считается напорным, ограниченным, фильтрационное течение стационарным, однофазным. В случае, когда значение числа Рейнольдса Re больше критического значения Re^, закон Дарси нарушается, и фильтрационное течение описывается двучленным нелинейным законом Форхгеймера grad h = -(А + Bv)v, где А=1/К, К -коэффициент фильтрации, v - модуль вектора скорости фильтрации v, В -константа пористой среды. Для вычисления чисел Рейнольдса используется формула Ф.И. Котяхова, Г.Ф. Требина Re(v) = 4-V2K-v-p/(m1,5 • jll), где k = K-|li/(p-g), jli - коэффициент вязкости, р - плотность, g - ускорение свободного падения, m - пористость.

В п. 2.2 строится итерационный алгоритм решения поставленной задачи без декомпозиции области. На i-ой итерации решается задача фильтрации, где область применения нелинейного закона фильтрации берется с предыдущего шага (G1"1). После решения область G'k определяется теми конечными элементами, на которых числа Рейнольдса Re (v) больше критического значения Re^. Итерационный процесс останавливается в случае выполнения условия G'"1 = Gk и достижения заданной точности.

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

В п. 2.4 приводятся результаты численных экспериментов. Рассматривался пятислойный пласт с различными толщинами слоев и коэффициентами фильтрации. Пласт считался непроницаемым за исключением кровли, на которой задавалось граничное условие 2-го рода, и противоположных участков боковой поверхности пятого слоя с граничными условиями 1-го рода. Коэффициент В выбирался таким образом, чтобы на границе перехода от линейного закона к нелинейному скорости, вычисленные с использованием законов Дарси и Форхгеймера, были близки. Результаты численных экспериментов приводятся в случаях, когда эти скорости отличались не более, чем на Av(%)=l%, 0.1%, 0.01%. Алгоритм апробирован на многопроцессорной вычислительной системе МВС-1000 при решении задач стационарной фильтрации в трехмерном неоднородном пласте с различным числом процессоров и гидродинамически несовершенных скважин.

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

В п. 3.1 дается постановка задачи. Задача решается в области G с D. Область D представляет собой пласт, ограниченный кровлей, подошвой, боковыми поверхностями и поверхностями интервалов вскрытия скважин V^. На участках Г, и Г2 границы области D задаются граничные условия первого и второго рода. Область решения G определяется в процессе решения с учетом депрессионных поверхностей в прискважинных зонах. Участки граничных условий первого и второго рода ГС1, ГС2 определяются положением верхней границы Гнб области G, при этом Гс с Г02, Г11б = Гн + Гб, где h > z на Гн Hh = z наГ6.

В п. 3.2 строится итерационный метод для решения задачи напорно-безнапорной фильтрации жидкости, основанный на методе декомпозиции области аналогичном, изложенному в п. 1.4. Для получения начальных значений решается задача на грубой сетке в предположении, что пласт является напорным. В каждой прискважинной подобласти решение представляется в виде суммы двух решений. Одно решение определено на сгущающейся сетке, другое - на дополнительной грубой сетке. Независимо решаются системы уравнений для сгущающихся участков сетки, затем совместно решается система уравнений для основной и дополнительных грубых сеток. При решении уравнений для сгущающихся участков сетки в прискважинных областях определяются участки верхней поверхности, на которой h < z, уточняются координаты z депрессионной поверхности из условия z = h, и все сетки перестраиваются от депрессионной поверхности. На дополнительных грубых сетках режим фильтрации считается напорным.

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

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

В п. 4.1 описывается схема метода решения уравнения насыщенности, основанного на декомпозиции области. Процесс построения алгоритма для определения насыщенностей на n-ом временном шаге состоит из следующих этапов:

1) Вычисляются по явной схеме фазовые расходы, выходящие из ячеек грубой сетки, окружающих ячейки сгущающихся участков прискважинных зон.

2) Определяются насыщенности по неявной схеме для ячеек сгущающихся участков, при этом фазовые расходы, входящие в прискважинные зоны, являются граничными условиями при решении указанных систем.

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

4) Подправляются насыщенности для ячеек грубой сетки.

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

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

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

1) использовался только метод с декомпозицией области по давлению;

2) использовался только метод с декомпозицией области по насыщенности;

3) использовались оба метода декомпозиции области.

Алгоритмы тестировались при решении модельных задач двухфазной фильтрации несжимаемой жидкости с различным числом скважин со сгущающимися сетками в прискважинных зонах на многопроцессорной вычислительной системе МКВС-Е112.

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

В п. 5.1 даётся постановка задачи. Рассматривается трехфазная изотермическая фильтрация нефти, воды и газа, подчиняющаяся линейному закону Дарси. Считается, что нефть и вода несжимаемы, и отсутствует массообмен между нефтяной и газовой фазами. Гравитационные и капиллярные силы не учитываются. В начальный момент времени считаются известными распределения давления и насыщенностей в пласте. На внешней поверхности пласта задаются граничные условия 1-го или 2-го рода для давления и известные доли фазовых компонентов на участках внешней поверхности, через которые поступают флюиды. На скважинах задаются забойные давления.

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

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

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

В заключении приводятся основные результаты и делаются выводы.

Апробация.

Основные результаты диссертационной работы докладывались и обсуждались:

-на XI международной конференции по вычислительной механике и современным прикладным программным системам (ВМСППС'2001, Москва-Истра, 2001);

- на IV научно-практической конференции молодых ученых и специалистов Республики Татарстан (Казань, 2001);

- на VIII Четаевской международной конференции « Аналитическая механика, устойчивость и управление движением» (Казань, 2002); на конференции «Современные проблемы гидрогеологии и гидромеханики» (Санкт-Петербург, 2002);

- на Всероссийской конференции «Высокопроизводительные вычисления и технологии (ВВТ-2003)»(Ижевск, 2003);

- на XVII сессии Международной школы по моделям механики сплошной среды (Казань, 2004);

- на международном семинаре «Супервычисления и математическое моделирование» (Саров, 2006);

- на итоговых научных конференциях КазНЦ РАН;

- на семинарах Института механики и машиностроения КазНЦ РАН. Основное содержание диссертации изложено в работах [11, 12, 22-30, 41].

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

Автор выражает глубокую благодарность своему научному руководителю канд.-физ.мат. наук Петру Алексеевичу Мазурову за предложенную тему исследований, постоянное внимание и поддержку при выполнении работы, а также сотрудникам группы математического моделирования Габидуллиной А.Н., Елесину А.В., Кадыровой А.Ш.

Работа выполнена в рамках программы фундаментальных исследований Президиума РАН "Параллельные вычисления на многопроцессорных вычислительных системах".

Положения, выносимые на защиту:

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

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

3. Численное тестирование алгоритмов с использованием различного числа процессоров и сгущающихся участков сетки.

5.4. Выводы.

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

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

Заключение.

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

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

Алгоритмы тестировались на многопроцессорных вычислительных системах МВС-1000 и МКВС-Е112 с различным числом процессоров и с1ущающихся участков сетки.

Показана эффективность построенных алгоритмов с декомпозицией области в сравнении с алгоритмами без декомпозиции.

1. Агошков В.И. Операторы Пуанкаре - Стеклова и методы разделения области в вариационных задачах / В.И. Агошков, В.И. Лебедев // Вычислительные процессы и системы. Вып. 2. - М.: Наука, 1985. - С. 173-227.

2. Азис X., Сеттари Э. Математическое моделирование пластовых систем / X. Азис, Э. Сеттари. М.: Недра, 1982, - 407 с.

3. Баренблатт Г.И. Движение жидкостей и газов в природных пластах / Г.И. Баренблатт, В.М. Ентов, В.М. Рыжик. М.: Недра, 1984, - 211 с.

4. Басниев К.С. Подземная гидравлика / К.С. Басниев, A.M. Власов, И.М. Кочина, В.М. Максимов. М.: Недра, 1986, - 303 с.

5. Бахвалов Н.С. Численные методы. Т.1. / Н.С. Бахвалов. -М.: Наука, 1975. -720 с.

6. Бахвалов Н.С. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. М.: Наука, 1987, - 632 с.

7. Вабищевич П.Н. Итерационные методы декомпозиции областей с налеганием для эллиптических краевых задач / П.Н. Вабищевич // Дифф. Уравнения. 1996. 32, №1. С. 923-927.

8. Воеводин В.В. Параллельные вычисления / В.В. Воеводин, Вл.В Воеводин.-СПб.: БХВ-Петербург, 2002. 608с.

9. Годунов С.К. Разностные схемы / С.К. Годунов, B.C. Рябенький. М.: Наука,1977. 440с.

10. Губайдуллин Д.А. Алгоритмы распараллеливания на сгущающихся сетках в задачах трехфазной фильтрации жидкости / Д.А. Губайдуллин, А.И. Никифоров, А.В. Цепаев // Вычислительные методы и программирование. -2007. Т.8. - №2. - С. 360 - 366.

11. Дмитриенко М.У. Оганесян JI.A. Вариант метода Шварца для прилегающих сеточных областей // Вычисления с разряженными матрицами. Новосибирск, 1981. С. 36-44.

12. Зенкевич О. Конечные элементы и аппроксимация / О. Зенкевич, К. Морган. М.: Мир, 1986,-318 с.

13. Ильин В.П. Методы неполной факторизации для решения алгебраических систем / В.П. Ильин.- М.: Физматлит, 1995. 288 с.

14. Кацнельсон В.Э., Меньшиков В.В. Об одном аналоге альтернирующего метода Шварца / В.Э. Кацнельсон, В.В. Меньшиков // Теория функций, функциональный анализ и их приложения.- Харьков: ХГУ, 1973.-Вып. 17.-С.206-215.

15. Коллинз Р. Течение жидкостей через пористые материалы / Р. Коллинз. М.: Мир, 1964.-350 с.

16. Копысов С.П. Методы декомпозиции и параллельные схемы метода конечных элементов / С.П. Копысов // Препринт ИПМ УрО РАН. Ижевск, 1999.

17. Кузнецов Ю.А. Вычислительные методы в подпространствах / Ю.А. Кузнецов. // Вычислительные процессы и системы М.: Наука, 1985. - С.265-350.

18. Лебедев В.И. Метод композиции / В.И. Лебедев. М.: ОВМ АН СССР, 1986.

19. Лебедев В.И., Агошков В.И. Обобщенный алгоритм Шварца с переменными параметрами / В.И. Лебедев, В.И. Агошков // Препринт ОВМ АН СССР. -М., 1981.

20. Мазуров П.А. Алгоритмы для распараллеливания решения задач двухфазной фильтрации жидкости на сетках со сгущающимися участками / П.А.

21. Мазуров, А.В. Цепаев // Тезисы международного семинара «Супервычисления и математическое моделирование». Саров. 2006. - С. 56.

22. Мазуров П.А. Алгоритмы для распараллеливания решения задач двухфазной фильтрации жидкости на сетках со сгущающимися участками / П.А. Мазуров, А.В. Цепаев // Вычислительные методы и программирование. -2006. Т.7. - №2. - С. 115-123.

23. Мазуров П.А. К решению задач фильтрации несжимаемой жидкости в трехмерных пластах с гидродинамически несовершенными скважинами / П.А. Мазуров, А.В. Цепаев // Математическое моделирование. 2002. Т. 14. №9. -С 121-123.

24. Мазуров П.А. Метод решения нелинейных задач фильтрации жидкости в трехмерных пластах с гидродинамически несовершенными скважинами / П.А. Мазуров, А.В. Цепаев // Математическое моделирование. 2004. Т. 16, №3, - С. 33-42.

25. Мазуров П.А. Решение трехмерных задач фильтрации жидкости на МВС-1000/16 на сетках со сгущающимися участками / П.А. Мазуров, А.В. Цепаев // Актуальные проблемы механики сплошной среды. Казань. 2004. - С.45-56.

26. Марчук Г.Е. Методы вычислительной математики / Г.Е. Марчук. М.: Наука, 1989. - 608 с.

27. Матеева Э.И. О разделении областей при решении краевых задач для уравнения Пуассона в областях сложной формы / Э.И. Матеева, Б.В. Пальцев // ЖВМ и МФ. 1973.- Т. 13, №6, - С.1441-1452.

28. Мацокин A.M. Связь метода окаймления с методом фиктивных компонент и' методом альтернирования по подпространствам / A.M. Мацокин // Диф.уравнения с частными производными. Новосибирск: Наука, 1986.

29. Мацокин A.M., Непомнящих С.В. Метод альтернирования Шварца в подпространствах. / A.M. Мацокин, С.В. Непомнящих // 1985. №10. С.61-66.

30. Мироненко В.А. Динамика подземных вод / В.А. Мироненко. М.: Изд-во Московского государственного горного университета, 1996. - 519 с.

31. Михлин С.Г. Об алгоритме Шварца / С.Г. Михлин // ДАН СССР. 1951. -Т.77, №4. - С.569-571.

32. Никольский Е.Н. Алгоритм Шварца в задачах теории упругости в напряжениях / Е.Н. Никольский // ДАН СССР. 1960. - Т.35, №3. - С.549-552.

33. Норри Д. Введение в метод конечных элементов / Д. Норри, Ж. де Фриз. -М.: Мир, 1981.-304 с.

34. Флетчер К. Численные методы на основе метода Галеркина / К. Флетчер. -М.: Мир, 1988.-352 с.

35. Цвик Л.Б. Обобщение алгоритма Шварца на случай областей без налегания / Л.Б. Цвик // ДАН СССР. 1975. - Е.224, №2.-309-312.

36. Bastian P. Numerical Computation of Multiphase Flows in Porous Media / P. Bastian. Habilitationsschrift, 1999. - 222 p.

37. Bramble J.H., Pasciar J.E., Xu J. Parallel multilevel preconditioners / J.H. Bramble, J.E. Pasciar, J. Xu // Math. Comput. 1990. - V.55. - P. 1-22.

38. Brenner S.C. BDDC and FETI-DP without matrices or vectors comput. Methods / S.C. Brenner, L.-Y. Sung // Appl. Mech. Engrg. 196. - 2007. - P. 1429-1435.

39. Cai X. Domain Decomposition in high-level parallelization of PDE codes / X. Cai // Proc. of the 11th Conf. on Domain Decomposition methods. Greenwich 1998. - P.388-395.

40. Cowsar L.C. Balancing domain decomposition for mixed finite elements / L.C. Cowsar, J. Mandel, M.F. Wheeler // Math. Comput. 1995. - V. 64. - P. 989-1015.

41. Dryja M. Towards a unified theory of domain decomposition algorithms for elliptic problems / M. Dryja, O.B. Widlund // Tech. Rep. Ultracomputer Note 167, Department of Computer Science, Courant Institute, 1989.

42. Farhat C. A scalable dual-primal domain decomposition method / C. Farhat, M. Lesoinne, K. Pierson // Numer. Linear Algebra Appl. 2000. - No.7. - P. 687-714.

43. Farhat C. FETI-DP: a dual-primal unified FETI method. I. A faster alternative to the two-level FETI method / C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, D. Rixen // Internat. J. Numer. Methods Engrg. 2001. - No. 50. - P. 1523-1544.

44. Fragakis Y. Force and displacement duality in Domain Decomposition Methods for Solid and Structural Mechanics / Y. Fragakis // Comput. Methods Appl. Mech. Engrg. -2007. No. 197. - P 584-599.

45. Gai Z. Multiplicative Schwarz algorithms for some nonsymmetric and indefinite problems / O. Widlund // SIAM J. Numer. Anal., 30. - 1993. - P. 936-952.

46. Huber R. Helmig R. Node centered finite volume discretizations for the numerical solution of multiphase flow in heterogeneous porous media // Computational Geosciences - 2002. - P. 141-164.

47. НШ M.C. Solving Groundwater Flow Problems by Conjugate-Gradient Methods and the Strongly Implicit Procedure / M.C. Hill // Water Resour. Res. -1990. -Vol.26. No.9. - P.1961-1969.

48. Irons В. M. A frontal solution program for finit-element analysis / В. M. Irons // Intern. J. on Numerical Methods in Engineering. 1970. - N. 2. -P.5-23.

49. Kopyssov S.P. Domain decomposition based on mechanical contact of substructures / S.P. Kopyssov, S.L. Ustyshanin // Intern. Conf. OFEA'2001 "Optimization of Finit Element Approximations & Splines Wavelets", St. Petersburg. 2001. -P.77-78.

50. Li B. Chen Z. Huan G. Control volume function approximation methods and their applications to modeling porous media flow // Advances in Water Resources. -2003.-P. 435-444.

51. Li J. FETI-DP, BDDC, and block Cholesky methods. / J. Li, O.B. Widlund // J. Numer. Methods Engrg. 2006. - No.66. - P.250-271.

52. Mandel J. Convergence of a balancing domain decomposition by constraints and energy minimization / J. Mandel, C.R. Dohrmann // Numer. Linear Algebra Appl.- 2003. No.10. - P.639-659.

53. Mandel J. Balancing domain decomposition / J. Mandel // Comm. Numer. Methods Engrg. 1993. - No.9. - P.233-241.

54. Mandel J. An algebraic theory for primal and dual substructuring methods by constraints /J. Mandel, C.R. Dohrmann, R. Tezaur // Appl. Numer. Math. 2005.- No.54. P.167-193.

55. Mazzia A. Putty M. Mixed-finite element and finite volume discretization for heavy brine simulations in groundwater // Journal of Computational and Applied Mathematics. 2002. - P. 191 - 213.

56. Novikov A.K. A parallel element-by-element conjugate gradient method with decreased communication cost / A.K. Novikov, S.P. Kopyssov // Intern. Summer School "Iterative Methods and Matrix Computations", Rostov-on Don. 2002. -P.450-454.

57. Przemieniecki J.S. Theory of matrix structural analysis. N/ Y.: VcGaw-Hill 1968.

58. Quarteroni A. Domain Decomposition Methods in Science and Engineering / A. Quarteroni, J. Periaux, Y. Kuznetsov, O. Widlund, eds. // Sixth International Symposium on Domain Decomposition Methods, June 15-19. 1992. - P.323-329.

59. Tamim M. Recent developments in numerical simulation techniques of thermal recoveiy processes / M. Tamim , J.H. Abou-Kassem. S.M Farouq AH. // Journal of Petrolium Science and Engineering. 2000. - P. 283-289.