«Метод конечных объемов для задачи конвекции-диффузии и моделей двухфазных течений ...»

Учреждение Российской академии наук

Институт вычислительной математики РАН

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

Никитин Кирилл Дмитриевич

Метод конечных объемов

для задачи конвекции-диффузии

и моделей двухфазных течений

05.13.18 – Математическое моделирование,

численные методы и комплексы программ

ДИССЕРТАЦИЯ

на соискание ученой степени

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

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

Василевский Юрий Викторович Москва – 2010 Содержание Введение................................... 4 Обзор используемой терминологии.................. Глава 1. Монотонная консервативная схема для задачи конвек­ ции-диффузии.............................. 1.1. Стационарное уравнение конвекции-диффузии.......... 1.2. Монотонная нелинейная схема конечных объемов на сетках с многогранными ячейками...................... 1.3. Сеточная система и свойства дискретного решения....... 1.4. Численные эксперименты...................... Глава 2. Численная модель двухфазной фильтрации в неодно­ родной пористой среде......................... 2.1. Уравнения двухфазной фильтрации................ 2.2. Схема, неявная по давлению, явная по насыщенности (IMPES) 2.3. Полностью неявная схема...................... 2.4. Схемы дискретизации потоков................... 2.5. Численные эксперименты...................... Глава 3. Численная модель течения несжимаемой жидкости со свободной поверхностью........................ 3.1. Математическая модель....................... 3.2. Численное интегрирование по времени.............. 3.3. Пространственная дискретизация на адаптивных сетках.... 3.4. Численные эксперименты...................... Заключение.................................. Литература.................................. Введение При решении современных инженерных и научных задач часто возни­ кает необходимость численного моделирования двухфазных течений. В на­ стоящей работе рассматриваются два подхода к представлению двухфазного течения: две фазы либо заполняют один объем, выражаясь при этом через насыщенности, либо имеют явную границу раздела.

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

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

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

Процесс решения задач математической физики можно разделить на несколько этапов: построение расчетной сетки, дискретизация [10, 11, 21], позволяющая преобразовать дифференциальное уравнение в систему алгеб­ раических уравнений, и решение системы алгебраических уравнений [9, 12, 24, 82]. В данной работе центральную часть занимают дискретизации.

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

В инженерном сообществе востребованы простые консервативные схемы, применимые на произвольных неструктурированных сетках для задач с ани­ зотропными свойствами среды. Кроме того, во многих практических задачах важно, чтобы численное решение отвечало определенным физическим требо­ ваниям, например, было неотрицательно. Консервативные линейные схемы на неструктурированных сетках хорошо известны: это метод конечных объ­ емов с многоточечной дискретизацией потока [27] (MPFA – multipoint flux approximaion), метод смешанных конечных элементов [39] (MFE – mixed finite element), метод опорных операторов [20, 61] (MFD – mimetic finite difference).

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

В первой главе диссертационной работы рассматривается метод конеч­ ных объемов для задачи конвекции-диффузии, который обеспечивает аппрок­ симацию потоков и сохраняет неотрицательность дискретного решения. Ме­ тод принадлежит к классу методов с нелинейной дискретизацией потока [1, 6, 46, 58, 59, 62–64, 72, 95]. Предлагаемый метод основан на схеме дис­ кретизации уравнения диффузии с полным анизотропным разрывным тензо­ ром диффузии на неструктурированных сетках с многогранными ячейками [4, 46]. В настоящей работе предлагается расширение данной схемы на случай уравнения конвекции-диффузии с разрывным полем скорости [72]. Стоит от­ метить, что вопрос эффективного решения систем уравнений, возникающих при дискретизации задачи конвекции-диффузии, не рассматривается в дан­ ной работе. Этому вопросу посвящен ряд публикаций [19, 30, 45, 50] Одной из основных трудностей при решении уравнений конвекции-диф­ фузии является подавление нефизичных осцилляций в численном решении.

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

В методах конечных элементов распространенным способом подавления нефизичных осцилляций является метод SUPG (Streamline Upwind / Petrov-Galerkin) [40]. Однако, осцилляции все равно могут появляться при реше­ нии задачи методом SUPG. Методы, уменьшающие нефизичных осцилляции, (SOLD – Spurious Oscillations at Layers Diminishing) [56] являются обобщением метода SUPG и удовлетворяют дискретному принципу максимума, по край­ ней мере, на некоторых модельных задачах.

Для дискретизации конвективных потоков можно использовать проти­ вопотоковую аппроксимацию, контролируемую через ограничение наклона [37, 42, 60] или внесение искусственной вязкости [35, 67]. Предлагаемая в диссертационной работе дискретизация конвективного потока является рас­ ширением предложенного в [64] двумерного метода на случай трехмерного пространства. Метод основан на идее монотонной противопотоковой схемы для законов сохранения (MUSCL – Monotone Upstream-centered Schemes for Conservation Laws) [93] и использует кусочно-линейное разрывное восполне­ ние с ограничителем наклона [54] для решения на многогранных ячейках.

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

Для дискретизации диффузионного потока применяется нелинейная двух­ точечная дискретизация потока на многогранных ячейках, предложенная в [46]. Идея монотонного метода конечных объемов для параболических уравне­ ний на треугольных сетках предложена К. ЛеПотье в [58]. В дальнейшем ме­ тод был распространен на более широкий класс сеток и уравнений [1, 6, 62, 95] и требовал интерполяции решения с основных переменных, определенных в ячейках, на вспомогательные переменные, определенные в узлах сетки. Ис­ пользование интерполяции влияет на точность решения, а также на скорость итерационного решения нелинейной системы. Была разработана двумерная схема, не зависящая от интерполяции [63], которая впоследствии была расши­ рена на трехмерный случай в [46]. Последняя схема формально не является безинтерполяционной и может потребовать интерполяции для небольшого числа вспомогательных переменных. Тем не менее, это не влияет на свойства схемы, поскольку бльшая часть интерполяций выполняется на основании физических принципов, таких как непрерывность полного потока на гранях сетки, по которым идет разрыв тензора диффузии.

Предлагаемый метод конечных объемов точен на линейных решениях.

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

Отметим, что для задачи диффузии с диагональным диффузионным тензо­ ром на кубических сетах шаблон сводится к традиционному семиточечному шаблону.

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

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

На основании предложенных схем дискретизации диффузионного и кон­ вективного потоков строится численная модель двухфазной фильтрации в пористой среде [5, 16, 23, 31, 43]. Для дискретизации по времени используют­ ся два наиболее популярных метода: метод, неявный по давлению, явный по насыщенности (IMPES-метод) и полностью неявный метод. Оба метода под­ разумевают использование дискретизации диффузионного потока Дарси на гранях ячеек. Стабилизация схемы осуществляется путем использования про­ тивопотоковой аппроксимации для значений насыщенности на гранях. Для сохранения второго порядка значения на гранях вычисляются на основании линейного восполнения с ограничителем наклона, аналогичного тому, кото­ рый разработан для конвективного потока в первой главе диссертации.

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

Цель второй главы диссертации – показать, что качество дискретизации диффузионного потока оказывает большое влияние на основные показатели моделирования разработки месторождения, такие как объем добычи, время прорыва и поведение водяного фронта. Производится сравнение двух схем дискретизации диффузионного потока с двухточечным шаблоном: традици­ онной линейной схемы и новой нелинейной. Многоточечная аппроксимация потока в данной работе не рассматривается в силу ее высокой арифметиче­ ской сложности и немонотонности [27].

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

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

Возможны разные подходы к эффективному и точному моделированию течений со свободной поверхностью. Они включают в себя методы, явно от­ слеживающие свободную границу [90, 91], и методы, основанные на неявном восстановлении поверхности [79, 88]. Методы конечных разностей [53], конеч­ ных объемов [51] и конечных элементов [33, 34] применяются для дискрети­ зации задачи по пространству.

Самые важные процессы, описывающие поведение жидкости, обычно происходят вблизи свободной поверхности, а сама поверхность может пре­ терпевать серьезные топологические изменения. В связи с этим распростра­ ненной практикой является использование адаптивных сеток, сгущающихся к текущему положению свободной границы [36, 51].

Большинство адаптивных технологий основаны на локально сгущающих­ ся треугольных (тетраэдральных) сетках и конечно-элементных дискретиза­ циях [36, 49], которые позволяют отслеживать сложные формы, возникаю­ щие при продвижении свободной поверхности. Недостатком подобного под­ хода для задач с постоянно меняющимся положением поверхности являют­ ся высокие вычислительные затраты, связанные с перестроением расчетной сетки, хранением и обработкой данных. Адаптивные декартовы сетки, напро­ тив, хорошо подходят для динамического сгущения или разгрубления сетки.

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

Дискретизации, основанные на сетках типа четверичное или восьмерич­ ное дерево активно используются при моделировании течений со свободной поверхностью [70, 80, 87], однако точные и эффективные схемы для подобных сеток все еще требуют дополнительного изучения.

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

Традиционная математическая модель, описывающая течение вязкой несжимаемой жидкости со свободной границей, основана на одновременном решении уравнений Навье-Стокса и уравнения функции уровня, описываю­ щего поведение свободной поверхности [85]. Помимо известных сложностей, связанных с расчетом течения жидкости, обработка свободной поверхности привносит дополнительные трудности. Во-первых, нахождение области, в ко­ торой решаются уравнения Навье-Стокса, само по себе должно быть частью вычислительного процесса. Во-вторых, свободная граница может претерпе­ вать серьезные изменения, такие как слияние и разделение фронтов. Для отслеживания подобных изменений требуется адаптивно сгущать расчетную сетку вблизи поверхности. В-третьих, для корректного воспроизведения сил поверхностного натяжения необходимо вычислять вектор нормали и локаль­ ную кривизну свободной поверхности. Это, в свою очередь, требует, чтобы функция уровня, описывающая положение свободной границы, определяла расстояние (со знаком) до поверхности.

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

Используются декартовы гексаэдральные сетки типа восьмеричное дерево, динамически сгущающиеся к поверхности в каждый момент времени. Для дискретизации операторов дивергенции, градиента и Лапласа используются компактные конечно-объемные и конечно-разностные схемы. Свойство рас­ стояния до поверхности восстанавливается путем решения дискретного урав­ нения Эйконала (так называемая операция реинициализации). Метод частиц [47] дополняет метод функции уровня и позволяет лучше сохранять объем жидкости.

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

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

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

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

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

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

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

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

Апробация работы. Результаты диссертационной работы докладыва­ лись автором и обсуждались на научных семинарах Института вычислитель­ ной математики РАН, Института прикладной математики РАН, Вычисли­ тельного центра РАН, Московского физико-технического института, Меха­ нико-математического факультета МГУ им. М. В. Ломоносова, Upstream Research Center of ExxonMobil corp. г.Хьюстон (США) и на следующих научных конференциях: конференция “Тихоновские чтения” (МГУ, Москва, ноябрь 2007); конференции “Ломоносов” (МГУ, Москва, апрель 2008, апрель 2010);

конференции “Ломоносовские чтения” (МГУ, Москва, апрель 2009, апрель 2010); конференция молодых ученых “Технологии высокопроизводительных вычислений и компьютерного моделирования” (СПбГУ ИТМО, С.-Петербург, апрель 2009); международная конференция “SIAM Geosciences 2009” (Лейп­ циг, Германия, июнь 2009); международная конференция “Computational Methods in Applied Mathematics: CMAM-4” (Познань, Польша, июнь 2010); меж­ дународные конференции “NUMGRID-2008” и “NUMGRID-2010” (ВЦ РАН, Москва, июнь 2008, октябрь 2010); международный научный семинар “Advances on Numerical Methods for Multiphase and Free Surface Flows” (ИВМ РАН, Москва, июнь 2009); международный научный семинар “Computational Mathematics and Applications” (Технологический университет Тампере, Тампере, Финляндия, сентябрь 2009).

Публикации. Основные материалы диссертации опубликованы в 7 пе­ чатных работах: 4 статьи – в рецензируемых журналах, входящих в перечень ВАК, [16, 17, 71, 72] и 3 статьи – в сборниках научных трудов [15, 18, 73].

Личный вклад автора. В совместных работах [18, 71, 73] вклад автора заключался в разработке вычислительного ядра технологии моделирования течения вязкой несжимаемой жидкости со свободной границей и проведении численных экспериментов.

В совместной работе [72] вклад автора заключался в разработке нели­ нейной схемы дискретизации уравнения конвекции-диффузии в трехмерном пространстве, программной реализации метода и проведении численных экс­ периментов.

Структура и объем диссертации. Диссертационная работа состоит из введения, обзора используемой терминологии, трех глав, заключения и списка литературы из 95 наименований. Диссертационная работа содержит рисунков и 12 таблиц. Общий объем диссертационной работы – 105 страниц.

Благодарности Автор диссертационной работы выражает глубокую признательность на­ учному руководителю Ю. В. Василевскому за продолжительную поддержку, ценные советы и плодотворное обсуждение вопросов. Автор благодарен С. Ю.

Малясову и В. Г. Дядечко из Upstream Research Center of ExxonMobil corp. за помощь в постановке задачи о практическом моделировании процесса двух­ фазной фильтрации в пористой среде. Автор также выражает благодарность М. А. Ольшанскому, И. В Капырину, А. А. Данилову, К. М. Терехову и многим другим за помощь в обсуждении идей и методов, используемых в диссерта­ ционной работе.

Работа над диссертацией была частично поддержана грантами РФФИ 08-01-00159-а, 09-01-00115-а, программой Президиума РАН “Фундаменталь­ ные науки – медицине”, федеральной целевой программой “Научные и научно­ педагогические кадры инновационной России”, а также грантом Upstream Research Center of ExxonMobil corp.

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

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

Каждая ячейка сетки является ячейкой звездного типа относительно центра масс, то есть каждая грань полностью видна из центра масс ячейки.

Аналогично каждая грань является плоской гранью звездного типа относи­ тельно центра масс грани.

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

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

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

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

Схему дискретизации уравнения конвекции диффузии, предлагаемую в первой главе будем называть монотонной. Под монотонностью метода или схемы здесь и далее подразумевается свойство сохранения неотрицательности получаемого численного решения при соответствующих граничных услови­ ях и правой части, обеспечиваемое монотонностью матрицы дискретизации дифференциального оператора [2, 22]. В общем случае дискретный принцип максимума [29] накладывает более строгие ограничения и требует сохране­ ния как минимума, так и максимума дискретного решения. Предлагаемый в данной работе метод дискретному принципу максимума не удовлетворяет.

Монотонная консервативная схема для задачи 1.1. Стационарное уравнение конвекции-диффузии Пусть – трехмерная многогранная область с границей, состоящей из двух частей: =, таких что =, и множество замкнуто и непусто: =, =. Предположим, что представима в виде объ­ единения конечного числа многогранных подобластей, = 1,...,, таких что = и = для =. На замыкании каждой подобласти определяется симметричный, положительно определенный, полный, возмож­ но анизотропный, непрерывный тензор диффузии K(x) с компонентами из 2 ( ) и непрерывное поле скорости v(x) с компонентами из 2 ( ), причем Рассмотрим стационарную задачу конвекции-диффузии для неизвестной концентрации : [7, 81]:

где 2 () – внешние источники, n – вектор внешней нормали, а и – заданные граничные условия. Обозначим через out – часть границы, на которой v · n 0, а через in = out. Предполагается, что out.

При вышеописанных условиях и соответствующих ограничениях на,, задача (1.1) имеет единственное обобщенное решение 2 () [3, 7].

Предполагаем, что выполняются достаточные условия неотрицательности ре­ шения () [3]: () 0, 0, 0. Физический смысл условий () и 0 состоит в отсутствии стоков внутри области и на границе.

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

1.2. Монотонная нелинейная схема конечных объемов на сетках с многогранными ячейками Для приближенного решения задачи (1.1) будет использоваться метод конечных объемов (также известный как метод интегральных тождеств [11] или метод интегрирования по подобластям [13]). Введем схему конечных объ­ емов с нелинейной двухточечной дискретизацией потока. Запишем уравнение баланса:

где q = K + v обозначает полный поток.

Пусть – конформная сетка, состоящая из многогранных ячеек с плоскими гранями. Предположим, что сетка связна по граням, то есть не мо­ жет быть разбита на две части, не имеющие ни одной общей грани. Обозначим через n внешнюю единичную нормаль к границе ячейки, а через n – вектор нормали к грани. Грани, принадлежащие границе расчетной обла­ сти, назовем внешними. Для них вектор n всегда внешний. Остальные грани назовем внутренними. Предположим, что |n | = | |, где | | обозначает пло­ щадь грани. Далее предположим, что сетка достаточно мелка и тензорная функция K(x) и поле скорости v(x) слабо меняются внутри каждой ячейки.

Тензор K и поле v могут иметь разрыв по граням ячейки, лежащим на, и менять направления (главные направления для K), однако потоки нормаль­ общей гранью.

Проинтегрируем уравнение (1.2) по многограннику и воспользуемся формулой Грина:

где q – средняя плотность потока на грани, а, – это либо 1, либо 1, в зависимости от взаимной ориентации векторов n и n.

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

где и – некоторые коэффициенты. В случае линейного метода данные коэффициенты фиксированы и равны между собой. Для нелинейного мето­ да они могут различаться и зависеть от значений концентрации в точках коллокации, окружающих + и. На грани представление потока совпадает с (1.4) с той лишь разницей, что одна из концентраций задана яв­ но. Для задачи Дирихле ( = ), подставив (1.4) в (1.3), получим систему из уравнений с неизвестными.

Таким образом, для метода конечных объемов со степенями свободы, ас­ социированными с ячейками, ключевым моментом является определение дис­ кретного потока (1.4). Предлагаемая схема дискретизации потока расширяет определение диффузионного потока [46] на случай конвективно-диффузион­ ного потока и основана на идеях, предложенных для двумерного случая в [64]. Для формирования дискретизации полного потока (1.4) разделим его на диффузионную и конвективную составляющие = + ±.

Введем обозначения. Пусть и – непересекающиеся множества внут­ ренних и внешних граней, соответственно. Подмножество внутренних гра­ ней включает грани, на которых происходит разрыв тензора диффузии K(x) или поля скорости v(x). Множество, в свою очередь, разбивается на под­ множества и, на элементах которых задаются граничные условия Дирихле и Неймана, соответственно. и обозначают подмножества внешних граней, принадлежащих границам out и in, соответственно. Нако­ нец, через и обозначаются множества граней и ребер многогранника, соответственно, а через – множество ребер грани.

Для каждой ячейки из точка коллокации x определяется в центре масс. Для каждой грани сопоставим точку коллокации с центром грани x. Кроме того, определим точки коллокации x в центрах Точки коллокации на гранях и ребрах являются вспомогательными. Зна­ чения концентрации в данных точках используются только в промежуточных вычислениях и не входят в итоговую систему алгебраических уравнений, хотя и могут влиять на значения коэффициентов в ней. Остальные точки колло­ кации будем называть основными. Значения концентрации в них образуют вектор неизвестных в итоговой системе уравнений.

Для каждой ячейки определим множество окружающих ее точек коллокации. Сначала включим в точку коллокации x. Затем для каж­ общую грань. Для оставшихся граней ( ) добавим точки коллокации x. Пусть ( ) обозначает количество элементов в множестве (мощность множества).

Аналогично, для каждой грани, принадлежащей ячейке, определим множество, окружающих точек коллокации. Инициализируем множество, = {x, x } и добавим к нему те точки из, которые соот­ ветствуют ячейкам или граням, у которых есть общие точки с. Мощность, обозначается через (, ).

Пусть для каждой пары ячейка-грань,, существуют три точки x,1, x,2 и x,3 из множества, для которых выполнено следую­ щее условие (см. Рис. 1.1): вектор конормали = K(x )n, отложенный из x, лежит внутри трехгранного угла, образованного векторами Коэффициенты,, определяются по следующим формулам:

где и |a b c| = |(a b) · c|.

Аналогично, предположим, что для каждой пары грань-ячейка,, : существуют три точки x,1, x,2 и x,3 из множе­ ства,, для которых вектор, = K (x )n, отложенный из x, лежит внутри трехгранного угла, образованного векторами и выполнены условия (1.6), (1.7). Здесь и далее под K (x ) будем понимать Тройки векторов t,1, t,1, t,1 будем для краткости называть триплета­ ми. Простой и эффективный алгоритм поиска триплета для пар и предложен в [46].

1.2.1. Нелинейная схема дискретизации диффузионного потока на внутренних гранях Для дискретизации диффузионного потока в трехмерном пространстве воспользуемся схемой, предложенной в [46].

Пусть – внутренняя грань, которую делят ячейки + и, и нормаль n внешняя для +, x± (или x± ) – точка коллокации в ячейке ±, а ± (или ± ) – значение дискретной концентрации в этой точке.

Рассмотрим случай и введем K = K(x ). Пусть = +. Исполь­ зуя введенные обозначения, определение производной по направлению и предположение (1.6), запишем Выведем схему дискретизации диффузионного потока с двухточечным шаблоном. Сперва заменим частные производные конечными разностями.

Затем выпишем другую аппроксимацию потока через грань, используя =. Введем индекс ±, чтобы различать q, для = + и =, а индекс опустим. Поскольку n является внутренней нормалью для, знак в правой части меняется:

где ±, ± и ± заданы формулой (1.7), и ±, обозначают концентрации в точках x±, из множества ±.

Определим новый дискретный поток как линейную комбинацию q · n с неотрицательными весами ± :

Веса выбираются таким образом, чтобы шаблон для q · n стал двухто­ чечным, и q · n аппроксимировало исходный диффузионный поток. Эти требования приводят к следующей системе:

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

Если точка коллокации x+, (или x, ), = 1, 2, 3, совпадает с точкой коллокации x (или x+ ), соответствующее слагаемое исключается из форму­ лы (1.13). Таким образом на регулярной кубической сетке шаблон принимает традиционный семиточечный вид: 6, 1, 1, 1, 1, 1, 1.

Решение (1.12) может быть выписано явно. Поскольку ±, ± и ± неот­ рицательны, при любых 0 получаем, что ± 0.

Если ± = 0, назначаем + = = 1/2. В противном случае, Подставив выражение для ± в (1.11) получим двухточечную формулу для диффузионного потока:

с неотрицательными коэффициентами Теперь рассмотрим случай, когда K+ (x ) и K (x ) различают­ ся. Рассмотрим концентрацию во вспомогательной точке коллокации x.

Получим двухточечные дискретизации потока отдельно для пар + и Неотрицательные коэффициенты +,,, получаются аналогично коэффициентам (1.15) на основании дискретных концентраций в точках кол­ локации из множеств ± и,±. Соответствующие конормальные векто­ ры к грани, ± = ±K± (x ) n являются внешними к ±. Непрерывность нормальных компонентов полного потока и поля скорости требует непрерыв­ ности нормального компонента диффузионного потока. Это предположение позволяет исключить из (1.16), (1.17) и вывести схему дискретизации диффузионного потока с двухточечным шаб­ лоном и коэффициентами Если оба коэффициента = 0, полагаем = ± /2 и = (+ + )/2.

Замечание 1.2.1. Отметим, что хотя формула (1.10) инвариантна от­ носительно изменения концентрации на постоянную величину, дискретный поток (1.14) определен корректно только для неотрицательных значений концентрации. Для анализа схемы в разделе 1.3 требуется расширить опре­ деление дискретного диффузионного потока на отрицательные концентра­ ции. Этого можно добиться путем добавления минимальной положитель­ ной константы ко всем значениям концентрации в (1.11), чтобы сделать их неотрицательными.

1.2.2. Нелинейная схема дискретизации конвективного потока на Предлагаемый метод построения дискретного конвективного потока был описан в [72] и является обобщением двумерного метода, предложенного в [64]. Для каждой внутренней грани конвективный поток дискретизируется с использованием линейного восстановления концентра­ ции на ячейке, взятой против потока где Определим восстановление как линейную функцию с вектором градиента g. Поскольку связана с центром масс ячейки, такое восстановление сохраняет среднее значение концентрации на ячейке при любом выборе g.

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

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

Затем множество получается из путем добавления дополнительных точек в случае, если точки множества лежат в одной плоскости. Если, отличные от x. Если = {x, x, x, x }, но объем тетраэдра об­ разованного этими четырьмя точками меньше | |, то также добавляем к элементы из, и, отличные от x.

Множество допустимых градиентов определяется через три ограни­ чения, предназначенные для подавления нефизичных осцилляций.

Во-первых, линейное восстановление, построенное на основании допусти­ мого градиента g, должно быть ограничено в точках коллокации x :

Из (1.23) следует, что g 0 для всех локальных минимумов и максимумов.

Во-вторых, в целях получения правильного знака конвективного потока, потребуем, чтобы восстанавливаемая функция была неотрицательна в точках Отметим, что когда центр грани x лежит вне выпуклой оболочки точек x, восстанавливаемая функция может принимать отрицательные значения в x даже при выполнении (1.23).

В-третьих, восстанавливаемая функция должна быть ограничена снизу во вспомогательных точках коллокации на out (они не принадлежат ):

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

Например, игнорирование ограничения (1.25) может порождать неустойчи­ вость при итерационном решении итоговой нелинейной системы алгебраиче­ ских уравнений.

Может быть доказано [64], что задача минимизации (1.22) с ограничени­ ями (1.23), (1.24), (1.25) имеет единственное решение.

Покажем, что из себя представляют ограничения, накладываемые на век­ тор градиента g R3. Ограничение (1.23) определяет полосу между двумя плоскостями в трехмерном пространстве:

Ограничения (1.24) и (1.25) определяют полупространства, отделенные плоскостями :

Обозначим через множество всех плоскостей { : x }, а через f – множество плоскостей.

Функционал отклонения в общем случае имеет изоповерхности эллип­ соидальной формы. Применим отображение (поворот и масштабирование), которое преобразует эллипсоиды в сферы (см. Рис. 1.2). Этот же опера­ тор отображает плоскости f в плоскости, точку g, являющуюся Рис. 1.2. Исходная задача с ограничениями и задача после преобразования.

решением задачи минимизации (1.22) без ограничений, в g. В результате задача минимизации с ограничениями сводится к проекции на пересечение отображений полос и полупространств. В алгоритме 1.1 для поиска решения g исходной задачи минимизации (1.22) с ограничениями используется реше­ ние преобразованной задачи.

Используя (1.20) и (1.21), получаем искомое представление конвективно­ го потока:

где ± состоят из линейной (первый порядок аппроксимации) и нелинейной (поправка второго порядка) частей:

индекс ± означает ±, и g± = g±.

Коэффициенты ± неотрицательны для положительной концентрации. Если = 0 в ячейке, тогда g должен быть нулевым вектором, и Алгоритм 1.1 Решение задачи минимизации (1.22) с ограничениями.

Методом наименьших квадратов найти решение g задачи (1.22) без ограничений;

если g удовлетворяет (1.23), (1.24) и (1.25) тогда Назначить g = {0, 0, 0}.

Применить отображение, чтобы преобразовать эллипсоидальные изо­ поверхности в сферические, и определить = (), g = (g ) и 10:

11:

12:

13:

14:

15:

16:

17:

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

18:

19:

20:

21:

22:

23:

24:

Применить обратное отображение g = (^ ).

25:

1.2.3. Дискретизация потоков на внешних гранях Рассмотрим внешнюю грань с граничным условием Неймана и ячейку, которой она принадлежит. Диффузионный поток через эту грань есть где, – среднее значение на грани. Можем считать грань ячейкой с нулевым объемом. Заменяя + и на и, соответственно, получим из формулы (1.26) схему дискретизации конвективного потока:

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

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

Схема дискретизации диффузионного потока представляется формулой где коэффициенты даны по формуле (1.15). Дискретизация конвективно­ го потока зависит от направления скорости на грани. Если, для дискретизации используются формулы (1.27) и (1.29). Если, исполь­ зуем где 1.2.4. Восстановление дискретного решения во вспомогательных Коэффициенты для диффузионного потока в (1.15), (1.19) и (1.28) могут зависеть от значений дискретного решения и во вспомогатель­ ных точках коллокации x, и x,. С другой стороны, дискретная система для метода конечных объемов формируется только для концентраций в основных точках коллокации. Значения,,, вычисляются по формулам (1.31), (1.32) в соответствии с граничными условиями Дирихле. Значения, восстанавливаются по формуле, должны восстанавливаться по имеющимся данным.

Значения концентрации на внешних гранях с граничным условием Ней­ мана восстанавливаются из концентраций по формулам (1.28) и (1.33).

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

1.3. Сеточная система и свойства дискретного решения Подставим двухточечную схему дискретизации (1.4) с неотрицательны­ ми коэффициентами = + ±, полученными из (1.15), (1.19) и (1.27), в уравнение баланса (1.3) и исключим концентрации на границе. Получим нелинейную систему из уравнений с неизвестными:

где – вектор дискретных концентраций в основных точках коллокации, а и – векторы значений концентрации во вспомогательных точках на гранях и ребрах, соответственно. Матрица M(,, ) составляется из матриц для внутренних граней и 1 1 матриц M (,, ) = (,, ) для внешних граней, на которых задано условие Дирихле. Вектор правой части G(,, ) строится по граничным данным и внешним источникам:

При () 0, 0 и 0, компоненты вектора G неотрицательны.

Для решения нелинейной системы (1.36) используется метод Пикара. Метод составления и решения алгебраической системы описан в Алгоритме 1.2.

Линейная система на Шаге 8 с несимметричной матрицей M(,, ) может быть решена, например, стабилизированным методом бисопряженных Алгоритм 1.2 Составление и решение нелинейной системы (1.36).

1: Для каждой пары ячейка-грань, и каждой пары грань­ ячейка, найти векторы t,1, t,2, t,3, удовлетворяющие условиям (1.5), (1.6) и (1.8), (1.6), соответственно.

ными значениями и порог non > 0.

Вычислить концентрации во вспомогательных точках коллокации на ребрах, используя (1.32) или арифметическое осреднение с граней, содер­ ) для каждой грани, используя (1.15), (1.19) и (1.27).

Вычислить вектор правой части G = G(,, ), используя (1.38).

Вычислить концентрации во вспомогательных точках коллока­ ции на гранях, используя (1.18), (1.28), (1.31), (1.33), и Вычислить концентрации во вспомогательных точках коллока­ 10:

ции на ребрах, используя (1.32) арифметическое осреднение с сосед­ 11:

градиентов с предобуславливателем (BiCGStab) [92]. Итерации метода оста­ навливаются, когда относительная норма невязки становится меньше lin.

Отметим, что, поскольку метод точен на линейных функциях, можно ожидать асимптотически второй порядок сходимости для концентрации и первый порядок для потоков.

Следующие две теоремы показывают, что решение (1.36) неотрицатель­ но при условии, что оно существует, и что значения концентрации на -й итерации метода Пикара неотрицательны, если lin = 0. Доказательства дву­ мерных аналогов теорем можно найти в [64].

и решение для (1.36) существует. Тогда 0.

Доказательство. Докажем от противного.

Рассмотрим ячейку с наименьшим значением концентрации и пред­ положим, что < 0. Без ограничения общности считаем, что векторы n внешние для, то есть = + в формулах для потоков. Поскольку ми­ нимально, концентрация на ячейке восполняется константой:. По определению конвективного потока:

где – соседняя ячейка, разделяющая грань. Напомним, что = + на внутренних гранях, и = на внешних гранях. Добавляя и вычитая, получим Из уравнения баланса (1.3) получим Имеем и, в предположении < 0, Из того, что минимально, следует (x ), а поскольку < 0, то, >. Пусть – вектор неотрицательных значений, полученный добавлением ко всем дискретным концентрациям положительной константы. Для получим Как отмечалось в Замечании 1.2.1, дискретный диффузионный поток для равен дискретному потоку для. Таким образом, q · n 0 и Поскольку 0, получим, что все члены в (1.39) неотрицательны и должны быть равны нулю. Следовательно что означает, что = для всех. Из этого следует, что вместо можно рассматривать любую соседнюю ячейку. Поскольку сетка связна по граням, делаем вывод, что - постоянна на. Из div v 0, следует, что =. Рассмотрим ячейку с гранью, и из (, ) = 0 получим, что неотрицательна. Это противоречит предположению и доказывает теорему.

и линейные системы в методе Пикара решаются точно, тогда 0 для Доказательство.

Покажем, что матрица M является М-матрицей при условии 0.

Значения концентраций во вспомогательных точках коллокации восстанав­ ливаются неотрицательными: 0 и 0. Кроме того, коэффициенты и ± по построению неотрицательны, следовательно = + ± 0.

Матрица M составляется суммированием матриц M (,, ) из (1.37) по всем граням сетки, следовательно ее диагональные элементы неотрицатель­ ны, внедиагональные – неположительны, сумма элементов в каждом столбце неотрицательна, а для граней с граничным условием типа Дирихле – сумма положительна. Поскольку расчетная сетка связна по граням, матрицы M и M неприводимы.

При описанных выше условиях в [94] доказывается, что матрица M является М-матрицей, и, следовательно, все элементы (M )1 положительны.

Поскольку транспонирование и обращение коммутируют, то все элементы (M )1 также положительны, и +1 = M1 G 0 при неотрицательной правой части G, что доказывает теорему.

Замечание 1.3.1. Теоремы 1.1 и 1.2 также справедливы и для линейной схемы дискретизации конвективного потока:

Эти утверждения позволяют нам называть предложенный метод моно­ тонным, хотя он и может не удовлетворять дискретному принципу максиму­ ма, см. раздел 1.4.4.

1.4. Численные эксперименты В экспериментах, посвященных анализу сходимости, полагаем =.

Для задач с доминирующей конвекцией это позволяет подбирать аналитиче­ ские решения такие, что все компоненты вектора правой части неотрицатель­ ны, G() 0, для любого неотрицательного вектора.

Будем использовать следующие дискретные 2 -нормы для измерения от­ носительных ошибок дискретизации для концентрации и полного потока q:

где | | – это среднее арифметическое объемов ячеек, содержащих грань.

Нелинейные итерации останавливаются, когда относительная невязка падает до non = 107. Критерием остановки для решения линейных систем является падение относительной невязки до lin = 1012. Для вычисления порядка сходимости будет использоваться алгоритм линейной регрессии.

Ниже рассмотрим три класса многогранных квазиравномерных сеток на единичном кубе [0, 1]3 (рис. 1.3).

Гексаэдральные сетки строятся на основе кубических сеток путем смеще­ ния узлов. Берется равномерная кубическая сетка с шагом. Узлы, лежащие в плоскостях = 0.5, = 0.5 и = 0.5, смещаются случайным образом в пределах окружности радиуса 0.2 на соответствующей плоскости. Осталь­ ные узлы сетки однозначно восстанавливаются по смещенным узлам, при условии, что грани ячеек плоские.

Призматические сетки имеют в основании квазиравномерную неструк­ турированную треугольную сетку и получаются тензорным произведением двумерной -сетки на одномерную -сетку. Обе сетки имеют одинаковый шаг. После построения начальной призматической сетки ее -плоскости на­ клоняются таким образом, чтобы они не пересекались, а высота ячеек была не меньше 0.75 и не больше 1.25.

Рис. 1.3. Примеры гексаэдральной (а), треугольной призматической (б), и тетраэдраль­ ной (в) сеток.

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

На рис. 1.3 представлены примеры сеток трех рассмотренных классов.

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

где I - единичная матрица.

Гексаэдральные Призматические Тетраэдральные Таблица 1.1. Ошибки сеточного решения и порядок сходимости для задачи с доминирую­ щей диффузией ( = 1).

Гексаэдральные Призматические Тетраэдральные Таблица 1.2. Ошибки сеточного решения и порядок сходимости для задачи с доминирую­ щей конвекцией ( = 0.01).

Источник и граничные условия Дирихле подбираются в соответствии с точным решением.

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

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

1.4.2. Анализ сходимости на задаче с разрывными тензором диффузии и полем скорости Далее рассмотрим сходимость к решению задачи с разрывными тензором диффузии и полем скорости. Пусть = (0, 1)3 разбита на две половины:

(1) = { < 0.5}, (2) = { > 0.5}, с границей, проходящей в плоскости = 0.5, с тензором K и скоростью v, претерпевающими разрыв по этой границе (см рис. 1.4).

и поле скорости v(x) = v() для x (), где Спектральное разложение K() = ( () ) () () демонстрирует боль­ шой скачок собственных чисел и направлений собственных векторов тензора K(x) в плоскости разрыва:

Рис. 1.4. Тензор K и скорость v, разрывающиеся вдоль плоскости = 0.5.

Рис. 1.5. Изолинии решения в плоскости для задачи с разрывными тензором диффузии и полем скорости.

Гексаэдральные Призматические Тетраэдральные Таблица 1.3. Ошибки сеточного решения и порядок сходимости для задачи с разрывными тензором диффузии и полем скорости.

Определим следующее точное решение задачи (1.1) с граничными условиями Дирихле = :

тогда правая часть равна Численные эксперименты проводились на гексаэдральных, призматических и тетраэдральных сетках, определенных выше. Сетки строились таким обра­ зом, чтобы граница = 0.5 точно приближалась гранями сетки. Изолинии решения в плоскости показаны на рис. 1.5. Результаты экспериментов, при­ веденные в таблице 1.3, показывают, что наличие разрыва в поле скорости и тензоре диффузии не влияет на порядок сходимости на всех рассматривае­ мых классах сеток.

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

[55, 56]). Следуя [56], возьмем Граничные условия Дирихле задаются следующим образом (см. рис. 1.6 (сле­ Для того, чтобы сохранить независимость решения от, введем однородные условия Неймана на границах, лежащих в плоскостях = 0 и = 1.

Рис. 1.6. Слева: схематичное представление условий задачи, справа: расположение обла­ стей 1, 2 и 3.

Точное решение имеет два пограничных слоя на границах, лежащих в плоскостях = 0 и = 1. Кроме того, имеется внутренний слой вдоль плос­ кости, проходящей по направлению поля скорости через линию (0, 0.7, ).

Численное решение неотрицательно во всех экспериментах, что соответ­ ствует Теоремам 1.1 и 1.2.

Вычисления проводятся на призматических сетках, в основании которых лежат структурированные (см. рис. 1.7а и 1.7б) и неструктурированные (см.

Рис. 1.7. Численное решение для сингулярно возмущенной задачи конвекции-диффузии ( = 1/32).

рис. 1.7в) треугольные сетки. Структурированные сетки получаются путем разбиения квадратных сеток северо-западной или северо-восточной диагона­ лями, соответственно. Эффективный шаг сетки составляет = 1/64 для экспериментов и = 1/32 для рисунков.

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

Зададим области 1, 2 и 3 (рис. 1.6 (справа)):

и 3 обозначает полосу ячеек, пересекающих линию = 0.25, Сначала введем величину (1.40), которая характеризует общую сумму осцилляций сеточного решения вне отрезка [0, 1] в области 1, покрывающей внутренний слой:

Далее введем величину (1.41), которая измеряет сумму осцилляций на пограничном слое, лежащем в 2 :

Наконец введем две величины (1.42) и (1.43), которые позволяют оценить толщину пограничного и внутреннего слоев, соответственно:

где Для непрерывного решения величины oscint = oscexp = 0, а smearint и smearexp зависят только от коэффициента диффузии, и для рассматриваемой задачи они должны быть очень малы. Для численного решения малые значения ве­ личин (1.40)–(1.43) характеризуют практически безосцилляционное и бездис­ сипативное дискретное решение. Чем меньше ширина внутреннего и погра­ ничного слоев в численном решении, тем меньший диссипативный эффект вносит рассматриваемая схема дискретизации.

Таблица 1.4. Сравнение осцилляций и диссипативности схем на сетке М1 (см. рис. 1.7а).

Таблица 1.5. Сравнение осцилляций и диссипативности схем на сетке М2 (см. рис. 1.7б).

Таблица 1.6. Сравнение осцилляций и диссипативности схем на сетке М3 (см. рис. 1.7в).

В таблицах 1.4, 1.5 и 1.6 приводятся результаты сравнения трех методов.

SUPG – традиционный конечно-элементный метод Петрова-Галеркина с про­ тивопотоковой стабилизацией – показывает малое размазывание на внутрен­ нем разрыве, однако вызывает заметные нефизичные осцилляции. Конечно­ элементный метод MH85, предложенный в [55], напротив, является практи­ чески безосцилляционным, однако вносит бльшую численную диффузию по сравнению с SUPG. По результатам сравнительного анализа, который прово­ дился в [56], метод MH85 был признан лучшим из более чем 20 методов по со­ вокупности введенных выше оценок. Предлагаемый в работе метод (FVMON) по интегральным оценкам, введенным в [56], не уступает лучшему участни­ ку данного сравнения. Бльшее значение ширины внутреннего разрыва на ‘северо-восточной’ призматической сетке вызвано бльшим шагом сетки в на­ правлении, перпендикулярном внутреннему слою.

Рис. 1.8. Изолинии решения в плоскости для задачи с условием непротекания на внеш­ ней границе.

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

Распространим тестовую задачу, исследованную в [27, 46], на уравнение конвекции-диффузии. Рассмотрим единичный куб с двумя вертикальными [3/11, 4/11] [5/11, 6/11], 2 = [7/11, 8/11] [5/11, 6/11]. Граница области разделяется на внешнюю часть, на которой задано однородное условие Неймана, и две внутренних части,1,,2 на поверхностях вырезов, где задаются граничные условия Дирихле: (x) = 0 для x,1 и (x) = для x,2.

Таблица 1.7. Максимальное значение концентрации для задачи с условием непротекания на внешней границе. = 1/22.

Задается анизотропный тензор диффузии где 1 = 3 = 1, 2 = 103, = 67.5, и () - матрица поворота в плоскости, ортогональной, на угол.

Поле скорости v = (102, 102, 0). Согласно дискретному принципу мак­ симума для эллиптических уравнений с частными производными, точное ре­ шение должно лежать между 0 и 1 и не может иметь строгих экстремумов на границе с условием непротекания.

Дискретное решение на кубической сетке с шагом = 1/22 показано на рис. 1.8. Оно неотрицательно, что отвечает теореме 1.1, но имеет экстремумы больше 1 вблизи границы области. Однако эти экстремумы быстро уменьша­ ются с измельчением сетки, как показано в таблице 1.7.

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

Численная модель двухфазной фильтрации в 2.1. Уравнения двухфазной фильтрации Рассмотрим уравнения двухфазной фильтрации в пористой среде [5, 31, 43, 77]. Фаза, которая смачивает среду больше, чем другая, именуется смачи­ вающей и обозначается индексом. Другая фаза называется несмачивающей и имеет индекс. В модели вода-нефть первая будет смачивающей, а послед­ няя – несмачивающей.

Запишем основные уравнения двухфазной фильтрации:

Закон сохранения массы в каждой фазе:

Закон Дарси:

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

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

В этих уравнениях – неизвестное давление фазы, – неизвестная насы­ щенность, u – неизвестная скорость Дарси, K – тензор абсолютной проница­ емости, – плотность, = () – вязкость, = () – фактор сжатия, = () – относительная проницаемость, = () – пористость, – гравитационный член, – глубина, а – источник или сток для скважины.

Граничные условия состоят из двух частей:

1. Условие непротекания (однородное условие Неймана) на границе рас­ 2. На скважинах задано фиксированное забойное давление bh.

Рассмотрим упрощенную модель скважин, в которой каждая скважина считается вертикальной с перфорацией в одном сеточном блоке. Источник или сток для скважины определяется по формуле [78]:

где – коэффициент продуктивности скважины, который не зависит от свойств жидкости, а определяется только свойствами среды.

В дискретных аналогах уравнений (2.1) и (2.2) для вычисления фазо­ давления на грани, а насыщенность аппроксимируются “вверх по потоку”:

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

2.2. Схема, неявная по давлению, явная по насыщенности (IMPES) Опишем метод интегрирования по времени системы уравнений двухфаз­ ной фильтрации, неявный по давлению и явный по насыщенности (IMPES – Implicit Pressure, Explicit Saturation) [5, 43]. Описываемая схема представляет из себя схему расщепления.

В качестве независимых переменных выберем давление нефти и насы­ щенность воды:

Определим полную скорость Дарси: u = u + u.

При фиксированной пористости среды и плотностях жидкостей имеем:

а значит:

Воспользовавшись (2.4), перепишем (2.2) в следующем виде:

где = + – полная мобильность.

Подставив (2.7) в (2.6), получим уравнение для давления:

Фазовые скорости u и u могут быть выражены через полную скорость u следующим образом:

Применяя (2.4) и (2.7) к уравнениям (2.1) и (2.2) (для = ), выведем уравнение для насыщенности:

Далее сформулируем метод IMPES как серию временных подшагов:

(1) Решить уравнение (2.8) и по текущему значению насыщенности получить текущее давление (неявная схема):

(2) Использовать (2.7) для нахождения текущей скорости Дарси u по текущим значениям насыщенности и давления :

(3) Явной схемой решить уравнение (2.9), чтобы по текущим насыщен­ ности, давлению и скорости получить насыщенность +1 на следу­ ющем временном шаге:

Заметим, что уравнение (2.10) представляет собой стационарное урав­ нение диффузии с тензором диффузии K и нелинейной правой частью.

При формировании диффузионных потоков в левой и правой частях уравне­ ния (2.10), а также в правых частях уравнений (2.11) и (2.12) воспользуемся линейной или нелинейной схемой дискретизации, которые будут введены в разделе 2.4. В разделе 2.5 приводится сравнительный анализ результатов ис­ пользования линейной и нелинейной схем дискретизации.

2.3. Полностью неявная схема Выведем полностью неявную схему для интегрирования по времени си­ стемы уравнений двухфазной фильтрации (2.1)-(2.4). Сначала воспользуемся неявной схемой для уравнений сохранения массы (2.1):

Теперь можно записать нелинейную невязку уравнений для -го прибли­ жения к величине, изменяемой на временном шаге + 1 в ячейке :

Дискретный аналог уравнения (2.13) может быть записан в виде:

для всех ячеек сетки в каждый момент времени.

Воспользуемся методом Ньютона для решения нелинейной системы (2.15) со скоростями Дарси (2.2):

где – номер шага метода Ньютона, – вектор главных (независимых) пере­ менных по всем ячейкам сетки:

– вектор нелинейных невязок по всем ячейкам сетки:

и – матрица якобиана:

Метод Ньютона останавливается, когда норма вектора невязок опускается ниже nwt.

Рассмотрим построение матрицы якобиана. Разделим невязки на две ча­ сти: аккумуляцию (включая скважины) и перенос,, =, +,, где Следует учитывать следующие зависимости параметров системы от глав­ ных переменных:

= ( ), где ( ) – заданная функция от насыщенности (см.

= ( ) – заданная функция от насыщенности, = ( ) – заданная функция от давления, = ( ) – заданная функция от давления, При формировании матрицы якобиана будем учитывать вклад вариации каждой из двух составляющих невязки.

2.3.1. Аккумуляция Рассмотрим вариацию аккумуляционной составляющей невязки:

Для производящей скважины имеем:

для нагнетательной скважины:

2.3.2. Перенос Теперь рассмотрим составляющую переноса в нелинейной невязке, осно­ ванную на потоках Дарси:

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

Здесь = (), – противопотоковая аппроксимация для значения на­ сыщенности на грани. С другой стороны, = () и = (), где Используя (2.20) и (2.21) получим следующее представление вариации потока для каждой из двух фаз:

где Считаем коэффициенты постоянными на каждой итерации метода Ньютона. В этом случае разница между линейной и нелинейной дискретиза­ циями потока заключается в способе вычисления для (2.22) и (2.23), но не в шаблоне для разреженной матрицы.

2.4. Схемы дискретизации потоков Для дискретизации потоков в (2.10)-(2.12) и (2.19) воспользуемся тради­ ционной линейной и новой нелинейной схемами дискретизации потоков Дар­ си.

2.4.1. Линейная схема Рассмотрим случай неортогональной сетки с анизотропным тензором проводимости: конормали K n и векторы t, соединяющие точки коллока­ ции, могут не быть ортогональными к граням сетки (см. рис. 2.1).

Считаем, что |n | = | | и обозначим ± = (x± ). Для потока через внутреннюю грань имеем:

Линейная двухточечная дискретизация t -компонента градиента давле­ ния имеет вид:

получим:

Потоки через внешние грани определяются через граничные условия Неймана.

В случае K-ортогональной сетки, когда K n и t сонаправлены, выраже­ ние (2.26) принимает вид центральной конечной разности и аппроксимирует поток как минимум с первым порядком точности. В общем случае линейная схема может не обеспечивать аппроксимации.

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

В случае K-ортогональной сетки, нелинейная схема принимает вид цен­ тральной конечной разности с постоянными коэффициентами и совпадает с линейной схемой.

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

2.5. Численные эксперименты Проведем сравнительный анализ линейной и нелинейной схем дискрети­ зации диффузионного потока. Рассмотрим псевдо-двумерные задачи на гек­ саэдральных сетках, состоящих 1 ячеек.

Здесь и далее величины приводятся в метрической системе мер. В скоб­ ках приводятся аналоги в британской системе. Приведем соотношения вели­ чин:

объем: 1 ббл (bbl, баррель) = 158.988 л. = 0.158988 м3, давление: 1 пси (psi, фунт на квадратный дюйм) = 6894.757 Паскаля, длина: 1 фт (ft, фут) = 0.3048 метра, проницаемость: 1 мд (md, миллидарси) = 9.869233 · 1016 м2, вязкость: 1 сп (cp, сантипуаз) = 0.01 г/(см · с) = 0.001 Н · с/м2.

Для экспериментов выбирались следующие свойства среды и жидкостей.

Относительные проницаемости показаны на рис. 2.2 (слева), зависимость капиллярного давления от водонасыщенности показана на рис. 2.2 (справа).

Вязкости и факторы сжатия заданы в таблице 2.1, и плотности равны = 9.797 · 103 Па/м (= 4.331 · 101 пси/фт) и = 8.817 · 103 Па/м (= 3.898 · 101 пси/фт). Фактор сжатия породы = 106.

4100 пcи), на производящей скважине – 2.689·104 кПа (= 3900 пcи). Индексы продуктивности скважины выбирались = 2.67 · 1012 м3 (= 10 день·пси ).

Относительная проницаемость Рис. 2.2. Слева: относительные проницаемости нефти и воды, справа: зависимость капил­ лярного давления от насыщенности.

2.5.1. Ортогональные сетки Численные эксперименты, представленные в данном разделе, показыва­ ют, что линейная и нелинейная схемы дискретизации потока дают одинако­ вые результаты на ортогональных сетках.

Рассматривается процесс заводнения в области 360.7 360.7 15.5 м, в двух углах которой расположены скважины: нагнетательная и произво­ дящая. Сетка, показанная на рис. 2.3, состоит из радиальной сетки вблизи скважин и неструктурированной (но все же K-ортогональной) сетки меж­ ду ними. Каждая скважина связана со всеми 48 угловыми ячейками с за­ данным коэффициентом продуктивности скважины = 3.05 · 1013 м (= 1.142 день·пси ). Среда изотропная. Тензор абсолютной проницаемости ска­ лярный K = I, где = 9.87 · 1014 м2 (= 100 мд).

Рис. 2.3. Неструктурированная ортогональная сетка. Водонасыщенность в момент времени = 1500 дней.

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

2.5.2. Неортогональные сетки Проведение сравнительного эксперимента на неортогональной сетке пред­ полагает следующий план действий: (1) построить регулярную сетку с кубическими ячейками; (2) подсоединить скважины и зафиксировать ячейки, к которым они подсоединены; (3) модифицировать оставшуюся часть сетки; (4) решить задачу с использованием линейной и нелинейной схем дис­ кретизации на модифицированной сетке; (5) сравнить результаты с теми, что были получены на регулярной сетке.

Рассматривается область 30.5 30.5 3.05 м, тензор абсолютной прони­ цаемости – скалярный, такой же, как в предыдущем эксперименте. Восполь­ зуемся следующей модификацией сетки: на равномерной сетке 64641 фик­ сируются первая и последняя строки, а так же первый и последний столбцы, с тем, чтобы угловые ячейки оставались без изменений. В центрах двух про­ тивоположных угловых ячеек помещаются скважины (рис. 2.4 слева). Далее вертикальная и горизонтальная центральные линии поворачиваются на 30 в направлении скважин (рис. 2.4 в центре) или на 30 от них (рис. 2.4 справа).

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

На рис. 2.5 показано распределение водонасыщенности на регулярной сетке Т1 (слева) и на сетке Т2 для линейной (в центре) и нелинейной (спра­ ва) схем дискретизации потока. На равномерной прямоугольной сетке линей­ ная и нелинейная схемы совпадают и имеют второй порядок аппроксимации потока. Расчет на регулярной сетке назовем контрольным. Хорошо видно, что при использовании линейной дискретизации, форма фронта отличается от формы контрольного фронта, в то время как для нелинейной дискрети­ зации отличия не заметны. Аналогично, на рис. 2.6 представлено сравнение результатов расчета на ортогональной сетке Т1 (слева) с расчетами на сет­ производящая скважина производящая скважина производящая скважина нагнетательная скважина нагнетательная скважина нагнетательная скважина Рис. 2.4. Пример ортогональной и неортогональных расчетных сеток.

Рис. 2.5. Водонасыщенность на сетках Т1 и Т2 для линейной и нелинейной схем в момент времени = 250 дней.

ке Т3 (в центре и справа). Как и в предыдущем случае, нелинейная схема дискретизации позволяет получить форму фронта, слабо отличающуюся от контрольной, в то время как для линейной схемы форма фронта не сохраня­ ется.

На рисунках 2.7 и 2.8 показаны динамики дебтов нефти и воды на про­ изводящей скважине. В момент, когда фронт обводнения доходит до про­ Рис. 2.6. Водонасыщенность на сетках Т1 и Т3 для линейной и нелинейной схем в момент времени = 250 дней.

Рис. 2.7. Динамика дебта нефти на производящей скважине на cетках Т1, Т2 и Т3.

изводящей скважины, производство нефти заметно падает. Данный момент называется временем водяного прорыва и является важным показателем мо­ делирования разработки месторождения. В табл. 2.2 показано время водяного прорыва при использовании линейной и нелинейной дискретизации на пред­ лагаемых расчетных сетках. Видно, что результаты, получаемые на неортого­ Рис. 2.8. Динамика дебта воды на производящей скважине на cетках Т1, Т2 и Т3.

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

2.5.3. Разрывный анизотропный тензор В следующем эксперименте будет рассматриваться регулярная сетка 64 1 из предыдущего примера, но тензор абсолютной проницаемости будет полным и анизотропным:

Расчетная область 30.5 30.5 3.05 м, 1 = 9.87 · 1013 м2 (= 1000 мд), 2 = 3 = 9.87 · 1015 м2 (= 10 мд), а углы поворота следующие (см. рис. 2.9):

нагнетательная скважина Рис. 2.9. Структура среды. Разрывный анизотропный тензор.

Стоит отметить, что при = 0 и = 90 линейная и нелинейная дис­ кретизации потока совпадают, различаясь лишь при = 45. Рисунок 2. иллюстрирует распределение водонасыщенности в момент времени = дней при использовании линейной (слева) и нелинейной (справа) дискретиза­ ций. На рисунке 2.11 показано поле давления в момент времени = 10 дней при использовании линейной (слева) и нелинейной (справа) дискретизаций.

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

Рис. 2.10. Водонасыщенность для линейной и нелинейной аппроксимаций потока в момент времени = 55 дней в задаче с разрывным анизотропным тензором.

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

2.5.4. Анализ вычислительной сложности Проанализируем влияние линейной и нелинейной схем дискретизации потока на сложность численного решения задачи двухфазной фильтрации на примере полностью неявного метода. Рассмотрим процесс заводнения на Рис. 2.11. Давление нефти для линейной и нелинейной аппроксимаций потока в момент времени = 10 дней в задаче с разрывным анизотропным тензором.

Рис. 2.12. Динамика дебтов нефти и воды для задачи с разрывным анизотропным тензо­ ром.

сетке Т3 из раздела 2.5.2 и сравним скорость работы метода Ньютона. Общее время расчета – 200 дней.

В таблице 2.3 показаны общее время расчета, общее число итераций ме­ Таблица 2.3. Шаг по времени (N – число шагов), время работы, число итераций, итераций на шаг, nwt = тода Ньютона и среднее число итераций на 1 шаг. Полностью неявный метод позволяет использовать большие шаги по времени без потери устойчивости.

Рассматриваются шаги от 1 до 20 дней. Критерий остановки метода Нью­ тона – nwt = 103. С ростом временного шага общее время расчета падает, однако число итераций, необходимых на каждый шаг, растет. При небольших шагах по времени общее время работы и число нелинейных итераций при ис­ пользовании нелинейной схемы дискретизации всего на 20% выше, чем для линейной схемы. С ростом шага по времени разница между схемами также возрастает (до 2.5 раз).

В таблице 2.4 указано общее время расчета, общее число итераций метода Ньютона и среднее число итераций на 1 шаг для другого критерия остановки метода Ньютона – nwt = 104. Видно, что при использовании нелинейной схемы невязка падает более медленно, и разница в числе итераций на данной задаче возрастает до 75% на малых шагах по времени и до 3.5 раз на больших.

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

Таблица 2.4. Шаг по времени (N – число шагов), время работы, число итераций, итераций на шаг, nwt = Выводы по второй главе Во второй главе рассмотрено применение новой нелинейной схемы дис­ кретизации диффузионного и конвективного потоков для модели двухфаз­ ной фильтрации. Результаты численных экспериментов показывают, что в случае неортогональных сеток и полных анизотропных тензоров абсолютной проницаемости нелинейная схема позволяет более точно и реалистично рас­ считывать распространение фронта обводнения и время водяного прорыва, чем традиционная двухточечная линейная схема.

Основные характеристики численной модели, созданной на основании но­ вой дискретизации, практически не зависят от неортогональности расчетной сетки и анизотропии тензора абсолютной проницаемости.

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

Численная модель течения несжимаемой жидкости со свободной поверхностью 3.1. Математическая модель Рассмотрим течение вязкой несжимаемой жидкости в области () R с границей, изменяющейся во времени. Предположим, что () = (), где – меняющаяся часть неподвижной границы, а () – свободная грани­ ца. На временном интервале (0, ] течение жидкости описывается уравнени­ ями Навье-Стокса для несжимаемой жидкости в области () где u - векторное поле скорости, - кинематическое давление, f - внешняя сила (например, сила тяжести), - плотность и - кинематическая вязкость.

В начальный момент времени = 0 расчетная область и поле скорости из­ вестны:

Предположим, что на неподвижной части границы поле скорости удо­ влетворяет граничному условию Дирихле:

где g задано. На свободной поверхности () накладывается кинематическое условие:

где n - вектор внешней нормали, а - нормальный компонент скорости на свободной границе (). Второе условие на () возникает при уравновешива­ нии нормального компонента тензора напряжения = [u + (u) ]/2 I силами поверхностного натяжения и внешним давлением ext :

Здесь – сумма главных кривизн, а – коэффициент поверхностного натя­ жения. Внешнее давление ext будем считать равным нулю, ext = 0.

Существующие подходы к численному решению (3.1)-(3.5) могут быть разделены на две группы: методы, основанные на отслеживании поверхности, и методы, в которых поверхность восстанавливается. Алгоритмы отслежива­ ния свободной поверхности основаны на явном представлении свободной по­ верхности, передвигающейся по закону (3.4). В работе используется алгоритм восстановления поверхности, основанный на неявном представлении () как нулевой изоповерхности глобально определенной функции (, x). Гладкая (по крайней мере, непрерывная по Липшицу) функция, для которой верно называется функцией уровня. Начальное условие (3.2) позволяет определить (0, x). Для > 0 функция уровня удовлетворяет следующему уравнению переноса [76]:

где u - любое гладкое поле скорости, для которого u = u на (). Использу­ емая математическая модель состоит из уравнений (3.1), (3.2), (3.3), (3.5) и (3.6). Отметим, что неявное определение () как нулевого уровня глобаль­ но определенной функции позволяет использовать численные алгоритмы, которые могут легко отслеживать сложные топологические изменения сво­ бодной поверхности. Функция уровня позволяет получать полезные геомет­ рические характеристики поверхности (). Например, единичная внешняя нормаль к () равна n = /||, а локальная кривизна поверхности рав­ на = · n. Для корректного вычисления кривизны и нормали к поверх­ ности необходимо, чтобы функция уровня являлась расстоянием (со знаком) до поверхности, то есть удовлетворяла уравнению эйконала 3.2. Численное интегрирование по времени Различные численные методы были предложены для интегрирования по времени уравнений (3.1),(3.6), начиная от полностью неявных схем и за­ канчивая методом дробных шагов. Полностью неявные схемы обеспечивают безусловную устойчивость, но имеют большую вычислительную сложность, поскольку требуют нескольких вложенных циклов итераций [52]. В данной работе применяется полунеявный метод расщепления [15, 71], который позво­ ляет избежать вложенных циклов, однако накладывает ограничения на шаг по времени. Алгоритм основан на хорошо известных процедурах расщепле­ ния, см. [44, 47, 75]. Стоит отметить, что ограничение на шаг по времени, накладываемое алгоритмом, на практике не является существенным.

Каждый временной шаг метода (при известных u(), (), () найти приближенные значения u( + ), ( + ), ( + )) состоит из подшагов, которые описывает алгоритм 3.1.

Остановимся подробнее на шагах алгоритма.

Алгоритм 3.1 Временной шаг метода.

1: Вычислить новое положение свободной поверхности.

Перенести функцию уровня по текущему полю скорости u().

Перенести частицы по текущему полю скорости u().

Выполнить коррекцию объема жидкости.

Восстановить расстояние (со знаком) до поверхности.

Перестроить сетку.

Вычислить новое поле скорости u( + ).

Перенести u( + ) по полю скорости u() с предыдущего шага.

Обновить u( + ) с учетом вязкости и внешних сил Спроектировать u( + ) на подпространство бездивергентных скоро­ 10:

Выбрать новый шаг по времени.

11:

3.2.1. Продвижение функции уровня: () ( + ).

На основании текущего поля скорости u() и текущего положения сво­ бодной поверхности () определим новое положение свободной поверхности ( + ).

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

2. Продолжим поле скоростей с поверхности на внешнюю часть объема жидкости: u()|() u()|R3. Из каждого узла сетки не принадлежа­ щего жидкости найдем ближайшую точку поверхности и возьмем зна­ чение скорости, проинтерполировав трилинейную функцию в данную точку поверхности. На практике продолжение необходимо выполнять лишь на небольшую часть внешнего объема: 0 < (x) < min.

3. Найдем ( + ) по формуле (3.6) численным интегрированием с ис­ пользованием полулагранжева метода [86] и продолженного поля ско­ рости. Полулагранжев метод состоит из следующих подшагов.

Сначала для каждой точки расчетной сетки y решим характеристиче­ ское уравнение назад по времени Характеристическое уравнение численно интегрируется со вторым по­ рядком точности:

Точка пространства x( + ) может не являться точкой расчетной сет­ ки, поэтому для вычисления u(x( + ), ) и u(x( + ), ) нужно использовать трилинейную интерполяцию для значения скорости.

Обозначим Для вычисления (x(), ) в (3.10) вновь используется интерполяция.

Обратим внимание, что на данном шаге свойство быть расстоянием (со знаком) до поверхности и сохранение объема могут нарушаться.

4. Применим коррекцию * ( + ) * ( + ) для того, чтобы выпол­ нялась глобальное сохранение массы;

5. Восстановим (реинициализируем) функцию уровня * ( + ) ( + ) так, чтобы ( + ) (приближенно) удовлетворяла уравнению эй­ конала (3.7).

Итоговая (+) неявно определяет новую расчетную область (+).

Для новой области обновим расчетную сетку.

3.2.2. Адаптивное перестроение сетки.

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

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

3.2.3. Динамика жидкости: {u(), ()} {u( + ), ( + )}.

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

1. Полулагранжев шаг выполняется, как было описано выше, с той лишь разницей, что y отвечает за узел, в котором задан определенный ком­ понент скорости, = 1, 2, 3, и уравнение (3.10) заменяется на Для вычисления (x(), ) в (3.11) вновь используется интерполяция для поля скорости.

2. Учет вязких и внешних сил:

3. Проекционный шаг: решая уравнение на давление ( + ):

обновляем скорость:

только на свободной поверхности, а следовательно, граничное условие в (3.12) определено корректно.

3.2.4. Адаптация шага по времени.

Выберем новый шаг по времени, удовлетворяющий условию Куранта:

где cfl - параметр, задающийся в приложении.

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

Замечание 3.2.1. Для более точного интегрирования (3.8) может быть по­ лезно разделить интервал времени [, + ] на частей и применить (3.9) на каждом подинтервале. При использовании = 2, 3 наблюдается заметное повышение точности вычислений по сравнению с = 1, однако возрастает и вычислительная сложность алгоритма.

3.3. Пространственная дискретизация на адаптивных 3.3.1. Расчетные сетки Для вычисления сил поверхностного натяжения и отслеживания слож­ ной топологии свободной поверхности, необходимо использовать расчетные сетки с мелким шагом в окрестности (). В этом случае использование рав­ номерных сеток становится непомерно дорогим, особенно в трехмерном про­ странстве. Локально измельчающиеся сетки требуют значительно меньших вычислительных затрат. Однако такие сетки должны динамически сгущаться и разгрубляться при продвижении свободной поверхности. При использова­ нии регулярных триангуляций (тетраэдризаций), перестроение сетки требует больших затрат оперативной памяти и большого числа переинтерполяций се­ точных данных. Данный шаг становится значительно менее затратным, если использовать декартовы сетки типа восьмеричное дерево с кубическими ячей­ ками. Двумерный аналог сетки типа восьмеричного дерева, сгущающейся к свободной поверхности, показан на рис. 3.1. Более подробную информацию о структуре данных для четверичного и восьмеричного дерева можно найти в [83, 84]. Использование кубических ячеек также полезно в связи с экономич­ ным распределением степеней свободы для скорости, а также эффективной и простой интерполяцией данных между двумя последовательными сетками, см. ниже. Кроме того, из-за наследственной иерархической структуры, для решения линейных систем могут применяться эффективные многоуровневые методы [32, 38].

Рис. 3.1. Слева: двумерное четверичное дерево, сгущающееся к свободной поверхности.

Справа: потеря точности отображения поверхности при переносе с мелкой сетки на более грубую.

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

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

рис. 3.1).

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

Для получения устойчивой аппроксимации мы используем разнесенное расположение неизвестных для скорости и давления при дискретизации урав­ нений Навье-Стокса [8, 53] (см рис. 3.2): давление аппроксимируется в цен­ трах ячеек, компоненты скорости - в центрах граней, функция уровня - в вершинах ячеек.

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

Рис. 3.2. Слева: Расположение степеней свободы на разнесенной сетке; - давление, {±, ±, ± } - компоненты скорости, - скалярная функция, например, функция уров­ ня. Справа: шаблон дискретизации (x ).

Во втором случае размеры ячеек, разделяющих грань, различные. Для простоты рассмотрим дискретизацию -компонента оператора градиента в центре грани x, как показано на рис. 3.2. Для этого рассмотрим центры четырех соседних ячеек x1,..., x4 и запишем формулу Тейлора для давления (x ) через значения (x ):

Пренебрегая членами второго порядка, мы получим систему из 4 линейных уравнений с 4 неизвестными.

Решая ее методом Крамера получаем следующий шаблон для -компонента оператора градиента:

где - определитель матрицы системы 4 4, a, - определители соответ­ ствующих миноров 3 3.

3.3.3. Оператор дивергенции скорости Для аппроксимации дивергенции скорости в центре x ячейки сетки воспользуемся формулой Гаусса где n - внешняя единичная нормаль к границе ячейки. Пусть - множе­ ство граней ячейки, т.е. =, и x обозначает центр. На основании (3.14) определим Благодаря разнесенному расположению точек, в которых хранятся компонен­ ты скорости, (u · n)(x ) определяется естественным образом.

3.3.4. Оператор Лапласа для скорости Определим u для заданной дискретной функции скорости u в уз­ лах для скорости (центрах граней) во внутренней части расчетной области.

Разнесенное расположение узлов для скорости используется для обеспечения устойчивости, однако это сильно осложняет построение компактного шабло­ на аппроксимации для оператора Лапласа на сетках типа восьмеричное дерево. По этой причине во многих работах, посвященных конечно-разност­ ным аппроксимациям уравнений Навье-Стокса на декартовых сетках типа восьмеричное дерево, см. [66, 70, 80], вязкие члены игнорируются (предел Эй­ лера), или степени свободы для скорости располагаются в узлах ячеек (это упрощает дискретизации, но требует введения дополнительной стабилизации, например, см. [74]). Для создания компактного шаблона дискретизации при использовании степеней свободы для скорости, расположенных в центрах гра­ ней, отдельно рассматриваются “нормальные” и “касательные” производные из для компонента скорости.

В качестве примера рассмотрим вычисление дискретного лапласиана для -компонента u. Разделим оператор Лапласа на две составляющие:

Сначала в центральной точке x каждой внутренней ячейки области, заполненной жидкостью, аналогично (3.14)–(3.15), вычислим первую произ­ водную компонента скорости по соответствующему направлению:

где n = (1, 2, 3 ) - внешняя нормаль к границе ячейки в точке x. Теперь значения 2 в -узлах внутренних ячеек находятся так же, как первый компонент градиента давления, то есть по формуле центральной разности, если грань разделяют две ячейки одного и того же размера, или по формуле (3.13) (где заменяется на ), если грань разделяют две ячейки разного размера.

Вклад 2 / 2 + 2 / 2 в формулу (3.16) вычисляется иначе. Рассмотрим любую -грань любой внутренней ячейки, заполненной жидкостью.

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

Если потоки q известны, то аналогично (3.15) считаем Остается вычислить q в центрах x ребер.

Для ячеек, окружающих ребро, может существовать три разных слу­ чая. (1) Ребро принадлежит двум -граням одинакового размера. (2) Рис. 3.3. Шаблон для потока через ребро в операторе Лапласа.