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

005534613

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

Терехов Кирилл Михайлович

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

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

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

1 О ОКТ 2013

Москва - 2013

/

005534613

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской

академии наук.

Научный руководитель:

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

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

Механико-математический факультет Федерального государственного бюджетного образовательного учреждения высшего профессионального образования "Московский государственный университет имени М.В.Ломоносова".

Защита состоится «31» октября 2013 г. в 15 часов на заседании диссертационного совета Д 002.045.01 при Федеральном государственном бюджетном учреждении науки Институт вычислительной математики Российской академии наук (ИВМ РАН), расположенном по адресу: 119333, г. Москва, ул. Губкина, д. 8, ауд. 127.

С диссертацией можно ознакомиться в библиотеке ИВМ РАН.

Автореферат разослан «ЗО сентября 2013 г. Ученый секретарь диссертационного совета,

доктор физико-математических наук Бочаров Г. А.

Ведущая организация:

Общая характеристика работы

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

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

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

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

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

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

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

Практическая значимость. Практическая значимость диссертацион-

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

На защиту выносятся следующие основные результаты:

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

2. С помощью данных алгоритмов разработан и реализован генератор сеток типа восьмеричное дерево.

3. На основе предложенного генератора разработана экономичная численная модель двухфазной фильтрации в пористой среде.

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

Апробация работы. Результаты диссертационной работы докладывались автором и обсуждались на научных семинарах Института вычислительной математики РАН, Института прикладной математики РАН, Вычислительного центра РАН, Института проблем безопасности развития атомной энергентики РАН, Upstream Research Center of ExxonMobil corp. (г.Хьюстон,

США) и на следующих научных конференциях: конференция "Тихоновские чтения", (МГУ, Москва, 2009 г.); конференция "Лобачевские чтения" (Казань, 2009 г.); 53-я научная конференция МФТИ (ИВМ РАН, 26 ноября 2010г.); международные конференции "Numerical geometry, grid generation and high performance computing" (ВЦ РАН, Москва, 13 октября 2010г., 17 декабря 2012г.); международная конференция "4th Workshop on Advanced Numerical Methods for Partial Differential Equation Analysis" (Санкт-Петербург, 22 августа 2011г.); международная конференция "Математическое моделирование природных катастроф и техногенных угроз" (Сьон, Швейцария, 20 августа 2013); европейская конференция "ENUMATH" (Лозанна, Швейцария, 26 августа 2013).

Публикации автора по теме диссертации. Основные материалы диссертации опубликованы в 10 печатных работах: 1 монография [1]; 5 статей - в рецензируемых журналах, входящих в перечень ВАК [5, 7-10]; 4 статьи -в сборниках научных трудов и материалов конференций [2-4, 6].

Личный вклад автора. В монографии [1] вклад автора заключался в предложении и реализации алгоритмов для хранения и работы с сетками общего вида, тестирования конкурентных пакетов; внедрения в программную платформу пакетов для решения систем линейных уравнений и пакетов для декомпозиции расчетной области на подобласти, приписанные к доступным процессорам; разработка и реализация программы для моделирования двухфазной фильтрации в пористой среде; параллелизация пакета "Povray" для визуализации посредством трассировки лучей. В совместной работе [9] вклад автора заключался в параллелизации существующей модели общей циркуляции океана, данный опыт лег в основу программной платформы для работы с сеточными данными. В совместных работах [5, 10] вклад заключался в разработке технологии моделирования течения вязкопластичной несжима-

емой жидкости со свободной границей, а именно в дискретизации оператора дивергенции от тензора напряжений, технологической цепочки для задания областей с реально!! топографией, верификации реализованного метода, постановке и проведении численных экспериментов со сходом оползня и разрушения дамбы. В совместной работе [8] вклад заключался в реализации динамических сеток тина восьмеричное дерево и полностью неявного метода для решения задачи двухфазной фильтрации в пористой среде. Был предложен критерий сгущения, разгрубления, метод интерполяции сеточных функций и произведена верификация метода. В совместной работе [7] вклад автора заключался в разработке конечно-разностного неявного метода для решения задачи конвекции-реакции-диффузии. Кроме того, в [7] автором была обнаружена неустойчивость, предложен и реализован метод стабилизации паразитного вихревого слоя; проведен ряд численных экспериментов для апробации метода и сравнения с референтными значениями.

Структура и объем диссертации. Диссертационная работа состоит из введения, обзора используемой терминологии, трех глав, заключения и списка литературы из 90 наименований. Диссертационная работа содержит 34 рисунка и 19 таблиц. Общий объем диссертационной работы - 124 страницы.

Содержание работы

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

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

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

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

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

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

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

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

• обмен данных и множеств элементов между процессорами,

• создание произвольного числа слоев фиктивных элементов,

• перераспределение элементов между процессорами,

а также множество промежуточных алгоритмов, позволяющих реализовать данные алгоритмы. Используемые алгоритмы в том числе описаны автором в монографии [1].

Во второй главе рассматривается численная модель, основанная на уравнениях Навье-Стокса. Предложены ппзкодиссипативные схемы дискретизации второго порядка точности по пространству и по времени для адаптивных сеток типа восьмеричное дерево с разнесенным расположением степеней свободы скорости. Данная схема является схемой расщепления, основанной на методе проекции Темама-Яненко-Шорина.

В разделе 2.1 приводятся уравнения течения вязкой несжимаемой жидкости:

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

— + и- Уи - г/Ди + Ур = 0 в П х (0,Т), от

(Цуи = 0 вПх[0,Г),

и|<=о = "о, р|«=о = Ро в О,

(1)

11|Г:

g

В разделе 2.2 предложена полунеявная схема расщепления для интегрирования уравнений (1) по времени:

аип+1 + рип + 7и1

п-1

-+

в О ((и" + £ (и" - и""1)) ■

= -Ур",

(2)

ип+1|г,

дп

= 0.

г2

Здесь £ = Дг"/Д^-\ а = 1+ £/(£ + 1). Р = + 1). 7 = + !)• Далее, спроектируем ип+1 на бездивергентное пространство, чтобы получить ип+1:

а(ип+1 - и"+1)/ДГ -Vq = 0, сИу и"+1 = 0,

п+1|

(3)

п-и"+1|Г1=п^, 9|г2 — 0-Задача (3) сводится к уравнению Пуассона для д>:

Д<7 = а(Иуип+1/А4п, = 0.

(4)

Г1

Получим новую скорость:

ип+1 = и„+1 + д^уд/а,

(5)

и новое давление:

„п+1 _ „тг

рп + д/а + ибху ип+1,

(6)

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

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

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

В разделе 2.5 предложен подход к заданию граничных условий на сетках типа восьмеричное дерево с учетом криволинейных границ. Данный подход

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

В разделе 2.6 приведены результаты следующих численных экспериментов, демонстрирующих корректность и эффективность предложенных дискретизаций:

• пример Эшера-Стейнмана с известным аналитическим решением для уравнений Навье-Стокса;

• моделирование течения в трехмерной каверне при разном числе Рей-нольдса (100,400,1000);

• обтекание цилиндров с квадратным и круглым сечением в узком канале при разном числе Рейнольдса (20,100,0-100).

На задаче Эшера-Стейнмана демонстрируется эффективность применения фильтра в подавлении паразитного вихревого слоя на сетках со сгущением и достижение схемой второго порядка точности. При моделировании течения в трехмерной каверне и обтекания цилиндров была продемонстрирована низкая численная вязкость схемы. Моделирование обтекания цилиндра с круглым сечением демонстрирует корректность подхода, описанного в разделе 2.5. Полученные данные, соответствующие моделированию течения жидкости в каверне, попадают в референтные значения, представленные в статье Вонга и Бэйкера. Коэффициенты подъемной силы и лобового сопротивления, полученные при моделировании обтекания цилиндров, сходятся к результатам, полученным Шафером и Туреком, Брааком и Рихтером.

Основные результаты, представленные в данной главе, опубликованы в статье [7].

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

В разделе 3.1 рассматривается течение двух фаз неперемешиваемых жидкостей в пористой среде. Фазу, которая обладает большей смачиваемостью, назовем смачивающей фазой и обозначим индексом т. Другую, несма-чивающую фазу обозначим индексом о. Здесь и далее 5„ и ра будут обозначать насыщенность и давление соответствующей фазы а = и>, о. Без потери общности будем называть фазу ги водой, а фазу о нефтью.

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

• Сохранение масс для каждой фазы:

д фБа = + а = 0 (7)

оЬ Ва

• Закон Дарси:

иа = —ЛаК (Vра - радУг), а = ги,о. (8)

Обе жидкости заполняют все пустоты:

5и, + 5,0=1- (9)

• Разница давлений между фазами определяется капиллярным давлением Рс = Рс{Зт)-

Ро Рт = Рс, (10)

Здесь К - абсолютный тензор проницаемости, ф - пористость, д - ускорение свободного падения, г - глубина, для фазы а: ра - неизвестное давление, Ба -неизвестная насыщенность, иа - неизвестные скорости Дарси, ра - неизвестная плотность, ра о - известная плотность фазы на поверхности, Ва = рао/ра - коэффициент объемного расширения, ца - вязкость, кта - относительная фазовая проницаемость, Ха = кта/(раВп) - мобильность, да - источник или сток.

Возьмем за основные неизвестные давление нефти р0 и насыщенность воды В дальнейшем мы также будем использовать следующие зависимости: кга = ¿Га(5ш), На = Ра(Ро), Ва = Ва(р0) И ф = ф0 (1 + СЕ(р0 - р°0)), где од - константа сжимаемости породы, ф0 и р°0 - некоторые заданные значения пористости и давления.

На границах резервуара поставим условия непротекания (однородные условия Неймана). Скважины учитываются через источники или стоки в уравнении (7). Предполагается, что каждая скважина является вертикальной, подключена к центру ячейки. Будем считать, что в скважине нет капиллярного давления. Таким образом, потоки нефти и воды в скважине зависят только от давления нефти. В данной работе используется формула для скважин, предложенная Писманом.

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

В разделе 3.3 описан конечно-объемный метод дискретизации по пространству, основанный на нелинейной двухточечной дискретизации дифузи-онного потока, предложенной Даниловым и Василевским.

В разделе 3.4 выписан метод вычисления матрицы частных производ-

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

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

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

В разделе 3.7 на задаче заводнения демонстрируется применение различных подходов учета вариации коэффициентов в нелинейной дискретизации дифузионного потока, предложенных в разделе 3-4-

В разделе 3.8 на данной численной модели демонстрируется применимость предложенных в первой главе алгоритмов для проведения параллельных расчетов на распределенных по процессорам сетках. Продемонстрирована эффективность распараллеливания модели на кластере ИВМ РАН и на суперкомпьютере ВМК МГУ В1иеСепе/Р.

Основные результаты опубликованы в статье [8].

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

Основные результаты и выводы

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

В работе предложена, реализована и проверена программная платформа для работы с сетками общего вида в параллельном режиме.

Предложена новая низкодиссинативная проекционная схема дискретизации уравнений Навье-Стокса для адаптивных сеток типа восьмеричное дерево с разнесенным расположением степеней свободы. Обнаружена неустойчивость в задаче проекции скорости на бездивергентное пространство и предложен низкочастотный фильтр для ее стабилизации. Предложенная схема показывает второй порядок сходимости в дискретных Ь2 и Ь^ нормах на аналитическом решении.

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

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

1. Ю. В. Василевский, И. Н. Конышш, Г. В. Копытов, К. М. Терехов. INMOST - Программная платформа и графическая среда для разработки параллельных численных моделей на сетках общего вида. Москва: Издательство Московского Университета, 2012. С. 144.

2. К. Д. Никитин, А. Ф. Сулейманов, К. М. Терехов. Технология моделирования течений со свободной поверхностью в реалистичных сценах // Труды Математического центра им. Н.И. Лобачевского. 2009. Т. 39. С. 305-307.

3. К. М. Терехов. Параллельная реализация модели общей циркуляции океана // Сборник тезисов лучших дипломных работ 2010. ВМИК МГУ, Москва: МАКС ПРЕСС, 2010. С. 30-31.

4. К. D. Nikitin, М. A. Olshanskii, К. М. Terekhov, Y. V. Vassilevski. Preserving distance property of level set function and simulation of free surface flows on adaptive grids // Численная геометрия, построение расчетных сеток и высокопроизводительные вычисления. 2010. Р. 25-32.

5. К. D. Nikitin, М. A. Olshanskii, К. М. Terekhov, Yu. V. Vassilevski. A numerical method for the simulation of free surface flows of viscoplastic fluid in 3D // Journal of Computational Mathematics. 2011. V. 29. P. 605-622.

6. K. D. Nikitin, M. A. Olshanskii, К. M. Terekhov, Yu. V. Vassilevski. Numerical modelling of viscoplastic free surface flows in complex 3D geometries // Proceedings of European Congress on Computational Methods in Applied Sciences and Engineering, ECCOMAS 2012. Vienna, Austria: September 10-12, 2012. P. 14. - 1 электрон, опт. диск (CD-ROM).

7. M. A. Olshanskii, К. M. Terekhov, Yu. V. Vassilevski. An octree-based solver

17

for the incompressible Navier-Stokes equations with enhanced stability and low dissipation. // Computers and Fluids. 2013. V. 84. P. 231-246.

8. K. M. Terekhov, Yu. V. Vassilevski. Two-phase water flooding simulations on dynamic adaptive octree grids with two-point nonlinear fluxes // Russian Journal of Numerical Analysis and Mathematical Modelling. 2013. V. 28, no. 3. P. 267-288.

9. K. M. Terekhov, E. M. Volodin, A. V. Gusev. Methods and efficiency estimation of parallel implementation of the sigma-model of general ocean circulation // Russian Journal of Numerical Analysis and Mathematical Modelling.

2011. V. 26, no. 2. P. 189-208.

10. Yu. V. Vassilevski, K. D. Nikitin, M. A. Olshanskii, K.M. Terekhov. CFD technology for 3D simulation of large-scale hydrodynamic events and disasters // Russian Journal of Numerical Analysis and Mathematical Modelling.

2012. V. 27, no. 4. P. 399-412.

Подписано в печать 27.09.13 г. Тираж 100 экз. Отпечатано в типографии «ГЕЛИОПРИНТ». Заказ № 341

Федеральное государственное бюджетное учреждение науки Институт вычислительной математики РАН

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

04201363309

Терехов Кирилл Михайлович

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

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

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель д. ф.-м. н., доцент Василевский Юрий Викторович

Москва - 2013

Содержание

Введение ........................................................................4

Обзор используемой терминологии ......................................18

Глава 1. Программная платформа для работы с сеточными данными ..........................................................................20

1.1. Операции модификации сетки........................................23

1.2. Адаптивные сетки типа восьмеричное дерево ......................24

1.3. Параллельные алгоритмы ............................................25

Глава 2. Численная модель течения вязкой несжимаемой жидкости ..............................................................................48

2.1. Математическая модель................................................48

2.2. Интегрирование по времени ..........................................49

2.3. Разложение Гельмгольца..............................................50

2.4. Дискретизация конвекции и диффузии..............................56

2.5. Расчетная область и граничные условия............................63

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

Глава 3. Численная модель двухфазной фильтрации в пористой

среде..........................................................................85

3.1. Математическая модель................................................85

3.2. Полностью неявная дискретизация..................................87

3.3. Конечно-объемный метод..............................................88

3.4. Метод вычисления Якобиана..........................................92

3.5. Сравнение линейной и нелинейной двухточечной аппроксимации

потока.................................. 97

3.6. Применение сеток типа восьмеричное дерево............101

3.7. Вычисление вариации нелинейной аппроксимации потока .... 106

3.8. Параллельный расчет.........................108

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

Литература...................................113

Введение

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

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

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

1 INMOST - Integrated Numerical Modeling Object-oriented Supercomputing Technologies

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

Детальный анализ подходов к хранению сеточной информации и иерархии связей между соседними элементами произвел Гаримелла в работе [71] Исходя из анализа, был выбран подход с оптимальным балансом между требуемой компьютерной памятью и сложностью вычисления всех необходимых связей

Существует ряд пакетов для работы с сетками, такие как MSTK2 [34], STK3 [28], МОАВ4 [82 83], FMDB5 [74] большинство из которых находится в разработке и по тем или иным причинам не удовлетворяют поставленным требованиям Некоторые из пакетов не предназначены для работы с динамически адаптируемыми сетками, некоторые пакеты предлагают недостаточный параллельный функционал, например, поддерживают всего один слой фиктивных элементов, некоторые в настоящий момент находятся в стадии активной разработки

Основной подход при параллельном решении уравнений математической физики является метод декомпозиции расчетной сетки и метод перекрытия сеток слоями фиктивных элементов [17] Метод декомпозиции расчетной сетки заключается в распределении исходной сетки между процессорами Задача распределения элементов сетки между процессорами оптимальным образом, для равномерной загрузки вычислительных узлов эквивалентна разрезанию связ-

2 MSTK - Mesh Toolkit сеточный инструментарий

3 STK - Sieira Toolkit

1 МОАВ - A Mesh-Oriented dat^Base сегочно-ориенхированная база данных

5 FMDB - Flexible distubuted Mesh DataBase шбкая распределенная сеючная база данных

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

Метод для параллельной работы с адаптивными сетками представлен в пакетах ЭТК и ЕМБВ. Он заключается в удалении слоя фиктивных элементов, перестроения сетки, а затем восстановления параллельного состояния новых сеточных элементов и слоев фиктивных элементов. Этот метод так же реализуем с помощью предложенных в этой работе алгоритмов, однако не является достаточно эффективным.

Во второй главе разработывается устойчивый низкодиссипативный метод решения уравнений Навье-Стокса, описывающих нестационарное течение вязкой несжимаемой жидкости. Сетки типа восьмеричное дерево завоевывают популярность в вычислительной механике и физике за счет своей простой прямоугольной структуры и вложенной иерархии. К примеру, такие динамически адаптируемые сетки были использованы в сочетании с конечно-объемными методами и методом Галеркина с применением к гиперболическим законам сохранения [31, 55, 70, 80]. Благодаря быстрому динамическому перестроению, такие сетки естественным образом подходят для моделирования задач с подвижными границами и течений со свободными поверхностями [33. 49, 50, 56, 69, 78].

Дискретизации для вязких и невязких уравнений течения жидкости уже были разработаны для динамических сеток типа восьмеричное дерево. Попи-нет [68] разработал конечно-объемную схему типа Годунова с использованием неразнесенных сеток, когда неизвестные компоненты скорости и давление расположены в центрах ячеек. Мин и Гибо [35, 52] разработали полу-лагранжев

метод, так же для неразнесенных сеток, но с неизвестными в вершинах сетки. В этих работах был применен специальный метод для стабилизации ложных мод давления, типичных для неразнесенных сеток. В работах [49, 50, 56] была использована MAC6 схема [3, 4, 40] с разнесенным расположением неизвестных, расширенная на сетки типа восьмеричное дерево.

Существует два преимущества разнесенного расположения неизвестных. Первое заключается в простом поэлементном наложении условия несжимаемости, выполнение которого эквивалентно сохранению массы. Второе заключается в устойчивости дискретизации по давлению, так как четные-нечетные ложные моды давления не могут быть представлены на такой сетке. Однако, такое расположение неизвестных усложняет построение схем высокого порядка, особенно в том случае, если рассматриваются разнесенные неизвестные на сетках типа восьмеричное дерево. К примеру, в работах [49, 50, 56] для конвекции был использован полу-лагранжев метод первого порядка.

В настоящей работе разрабатывается схема второго порядка точности, основанная на методе проекции Темама-Яненко-Шорина [5, 21, 36]. Для линеаризованной [63, 65] конвекции используются конечно-разностные противопоточ-ные схемы второго и третьего порядка с низкой численной вязкостью. Дискретизация основана на линейных и кубических интерполяциях по двум переменным, за счет чего шаблон дискретизации остается компактным, а матрица линейной системы при неявной дискретизации членов диффузии и конвекции остается разреженной. Для дискретизации задачи конвекции-диффузии по времени используется формула обратных разностей второго порядка [13]. За счет неявного шага конвекции удается избежать ограничения по Куранту на шаг по времени [22], так как это ограничение оказывается довольно сильным при расчете на

6 MAC - Marker and Cell

адаптивных сетках.

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

Одной из фундаментальных проблем сеток типа восьмеричное дерево является ступенчатая аппроксимация криволинейных границ. Различные подходы для аппроксимации краевых условий на криволинейных границах для прямоугольных сеток были предложены в работах [8, 30, 78, 86]. В этой работе предложен метод аппроксимации краевых условий типа Дирихле на криволинейных границах, применимый на сетках типа восьмеричное дерево.

Разработанный метод был проверен на ряде задач: аналитическое течение типа Бельтрами [29], течение в каверне с подвижной границей [43, 75, 90] и обтекание цилиндра в узком канале [19, 73].

Описанный метод является частью разрабатываемой модели, используемой для моделирования течений вязкопластичной жидкости со свободной поверхностью [60, 61]. Указанная модель была успешно применена к моделирова-

нию катастроф [89].

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

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

Существует несколько подходов к дискретизации уравнений двухфазной фильтрации воды и нефти по времени: IMPES7 [12, 76, 79] - условно устойчивый полунеявный метод и полностью неявный метод [26], обладающий безусловной устойчивостью. В этой работе используется полностью неявный метод, позволяющий избежать ограничения на шаг по времени при сгущении сетки.

Ранее Сухиновым [6] и Саадом [72] был предложен подход решения задачи фильтрации на адаптивных сетках типа четверичное дерево (в плоскости). Одним из недостатков подходов, предложенных в рассматриваемых работах, было

7 IMPES - IMplicit Pressure Explicit Saturation, неявное давление, явная насыщенность

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

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

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

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

В настоящей работе используется нелинейная конечно-объемная схема, суть которой заключается в использовании монотонной нелинейной двухточечной аппроксимации потока. Метод был впервые применен для параболических уравнений на треугольных сетках К. Лепотье [45]. Этот подход был применен

к широкому кругу задач [23, 47, 48, 58, 77, 88]. Метод позволяет работать с не К-ортогональными сетками и произвольными многогранными сетками.

Альтернативой нелинейной схеме является многоточечная схема аппроксимации потока [9]. Многоточечная схема является линейной, аппроксимирует концентрации со вторым порядком, но является только условно устойчивой [42] и условно монотонной [62].

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

Корректность и эффективность предложенных методов для задачи двухфазного заводнения демонстрируется на ряде численных экспериментов.

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

В первой главе рассмотрен подход к хранению сеток общего вида. На основе этого подхода описан алгоритм динамического измельчения сеток типа восьмеричное дерево.

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

В третьей главе рассмотрена задача двухфазной фильтрации в пористой

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

Актуальность работы. При численном моделировании задач математической физики часто приходится сталкиваться с недостатком компью�