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

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

Григорьев Александр Виссарионович

Численное моделирование фильтрации в трещиновато-пористых средах на основе модели двойной пористости

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

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

Якутск 2013

28 НОЯ 2013

005541102

Работа выполнена на кафедре прикладной математики и в Центре вычислительных технологий Института математики и информатики Северо-Восточного федерального университета имени М. К. Аммосова

Научный руководитель: Вабищевич Петр Николаевич, доктор

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

Защита состоится 20 декабря 2013 г. в 16ю часов на заседании диссертационного совета Д 212.306.04 при Северо-Восточном федеральном университете имени М. К. Аммосова, расположенного по адресу: 677000, г. Якутск, ул. Белинского 58, зал заседаний Ученого совета.

С диссертацией можно ознакомиться в научной библиотеке Северо-Восточного федерального университета имени М. К. Аммосова. Автореферат разослан 19 ноября 2013 года.

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

физико-математических наук, профессор, профессор Казанского Государственного Университета, г. Казань Мордовской Сергей Денисович, доктор технических наук, доцент, заведующий кафедрой ИТ Института Математики и Информатики СВФУ, г. Якутск

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

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

д.ф.-м.н.

Саввинова Н.А.

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

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

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

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

Математические модели движения жидкости в такой среде были разработаны в конце 50-х годов Г.И.Баренблаттом, Ю.П.Желтовым, И.Н.Кочиной. В современной литературе эта модель двойной пористости (в порах и в трещинах) известна как модель Баренблатта. Модель характеризуется наличием обмена давлениями между фазами.

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

Эрозия почвы — разрушение и снос верхних наиболее плодородных горизонтов почвы. Эффективная защита почв от водной эрозии возможна при плановом и систематическом внедрении комплекса противоэрозионных мероприятий, разработанного с учетом конкретных природно-географических условий каждого района или хозяйства. Также такие задачи ставятся при строительстве помещений, технических сооружений и транспортных путей сообщения. Защита естественных и техногенных массивов дисперсных грунтов от размыва является одной из самых актуальных проблем на территории России. На данный момент не существу-

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

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

• Построение вычислительных алгоритмов расщепления, реализующих модели двойной пористости;

• Создание программ для численного моделирования фильтрации в трещиновато-пористых средах;

• Тестирование разработанных вычислительных алгоритмов и программ на модельных тестовых двумерных задачах;

• Численное моделирование задачи просачивания жидкости в грунт при наличии в нем системы трещин.

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

• Построены и исследованы схемы расщепления для моделей двойной пористости;

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

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

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

• International Young Scientists Conference on Mathematical Modeling (Liniy, China, May 24-25, 2010);

• Supercomputer Technologies of Mathematical Modeling (Yakutsk, November 2830, 2011);

• Математическое моделирование развития северных территорий Российской Федерации (Якутск, 21-26 мая, 2012);

• 5th International Conference on Porous Media & Annual Meeting (Prague, Czech Republic, May 21-24, 2013);

• Supercomputer Technologies of Mathematical Modeling (Yakutsk, July 8-11, 2013).

Работа выполнена при поддержке грантов:

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

окружающую среду регионального конкурса ДАЛЬНИЙ ВОСТОК Российского фонда фундаментальных исследований;

• Разработка прикладного ПО для численного моделирования добычи углеводородного сырья на высокопроизводительных вычислительных системах №5542 ГЗ МО;

• грант 02.740.11.0041 Разработка симуляторов экологически безопасных технологий разработки и мониторинга месторождений ископаемых Арктики и регионов Севера Федеральной целевой программы.

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

Структура и объем диссертации. Диссертация состоит из введения, трех глав, заключения и списка литературы.

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

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

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

-д- — 7— П1У graa и — шу grad и — 0.

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

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

Рис. 1: Распределение давления в трещинах (сверху), в порах (снизу) при£ = 1. сти, отметить, что фильтрация в трещинах протекает быстрее чем в порах.

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

В задаче Коши для эволюционного уравнения первого порядка ищется функция и{Ь) € Н, удовлетворяющая уравнению

с[и

В— + Аи = /(4), г > О, по заданному /(£) е # и начальному условию

«(0) = и0.

Линейные операторы А и В — положительные, самосопряженные и стационарные. В случае псевдопараболической модели оператор

В = Е + -у А.

Для аддитивных разностных схем характерно разбиение оператора А на сумму операторов более простой структуры:

V

А = ^Аа, Аа = А*а> 0, а= 1,2,...,р.

а=1

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

{Е + 1Аа)~^- + Ааиа = Ш, 4 > 0, а = 1,2, ...,р

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

В работе строятся векторные аддитивные схемы. От исходного уравнения переходим к системе уравнений, когда компоненты вектора и = {щ, и.2, ...¡ир} определяются из

^+7 +ЕЛ»"/» = /(')• ¿>0,

0=1 0=1

иа(0) = и°, а =1,2,..., р. Построенная схема расщепления для векторной задачи имеет вид

(Е + 7Аа) (вУ*+1~У° + (1 - + £ (Е + \

V Т Т ' 3=1 Т

i п4-1 , п и , п—1 р „,п i

, Л Ур +Ур , Л Уа + 2Уа + Уд , V" Л УР + УР _ ,„п + -%- "-4- ^ Р-2- ^ '

0=1 Р=а+1

п = 0,1,..., а=1,2,...,р. Вычислительная реализация схемы связана с решением сеточных задач

(Е + (7б + I) Д>) = Х2, « = 1,2, ...,р.

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

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

- сИу(^1(а:)§гас1 «1) + г{х){их - и2) = !\{х,г),

ди

с2(х)-~ <Иу(<12(х) §гас1 и2) + г(х)(и2 - щ) = /2(ж, £).

где индексами 1 и 2 обозначены соответственно параметры трещиноватой и пористой сред, г(х)(и1 — и2) — обменный поток между средами, с[а — проницаемости, ¡1 — вязкость флюида, а = 1,2. Краевая задача рассматривается в ограниченной области П при 0 < t <Т.

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

с [ У1+ + / £гас! г/"+1 у^х + г / (у?+1 - = 0,

Jíl т Уп Ja

10

[ ^-— + й [ gгad у%+1 ^ у2йх + г [ (т/"+1 - = 0.

Jn т ¿п ¿п

Обезразмеривание проведено при постоянных коэффициентах и однородных правых частях так, что

С\ = С, (¡1 = 1, С2 = 1, ¿2 = <1.

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

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

с [ ~ + [ grad у"+1 grad у^х + г [ (у?+1 - У^хйх = 0,

Jn т Ju J п

[ ~ УК2с1х + й [ grad уТ1 gгad у2йх + г [ (у%+1 - у1+1)^х = 0. ./П Т Jn Jn

Для параллельной схемы расщепления имеем:

с [ У"+ ~ у7\)Хйх + [ gгad у"+1 gгad у^х + г [ (у?+1 - У2)у^х = 0,

[ У*+1 ~ У*у2йх + й [ grad у2+г ё^ас! у2йх + г [- = 0.

Уп т Jn ,/п

Обобщенная модель двойной пористости основывается на системе уравнений

Си{х)~- + с12(х)^- - div(dll(a:) grad щ) - Аш(с1п(х) grad и2) +г(ж)(и1 -и2) = Мх,Ь),

с2\{х)-^- + с22(х)-^- - (1^(^21 (ж) grad щ) - <Иу((122(х) gгad и2)

+г{х){и2 - щ) = /2(х,Ь). Построены схемы расщепления, которые основаны на явно-неявных аппроксимациях по времени. Рассмотрен, в частности, случай зацепления уравнений не только по младшим и старшим коэффициентам системы, но и зацепление по коэффициентам при производных по времени. Доказана безусловная устойчивость явно-неявных схем в соответствующих гильбертовых пространствах.

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

дв С ) д

1711 ЭрГ т ~ ^(^(Р!) 8гаЛ (Р1 + *з)) + г(х)(р! - р2) = /ь

т2 д^2 д* - ¿1у{К2{р2) gгad (р2 + Жз)) + г(х)(р 2 - л) = /2,

где индексами 1 и 2 обозначены параметры трещиноватой и пористой сред соответственно, ра = ра/рд — приведенное давление, ра — давление, та — пористость, ¿а(Ра) — насыщенность, Ка(иа) — гидравлическая проводимость, г(х)(щ — и2) — обменный поток между средами. Данные уравнения являются существенно нелинейными из-за нелинейных функций за(ра), Ка(ра). Вклад гравитационный силы в процесс просачивания воды осуществляется за счет последнего слагаемого в диффузионных операторах ёг\г(.ЙГа(ра) (ра + хл)).

На модельной двумерной задаче исследуется зависимость точности приближенного решения от физических и вычислительных параметров. На рис. 2 приводятся характерные результаты по распределению давления при шаге по времени г = 0.1 на расчетной сетке из 1132 треугольных элементов со сгущением в участках с наибольшими градиентами решения.

Рис. 2: Распределение насыщенности в трещинах (сверху), в порах (снизу) при

í = 20 ч.

Основное внимание уделяется численному моделированию задачи просачивания из канала с учетом системы колодцев. Трехмерные расчеты выполнены при использовании конечно-элементных аппроксимаций по пространству и шаге по времени г = 0.1 на расчетной сетке из 71574 тетраэдральных элементов. На рис. 3, 4 показано давление в трещинах и в порах.

Рис. 3: Распределение насыщенности в трещинах при I = 10 ч.

Рис. 4: Распределение насыщенности в порах при Ь = 10 ч.

Расчеты выполнены на вычислительном кластере СВФУ. На рис. 5 отображено разбиение области на не пересекающиеся подобласти. Автор применяет

Потоки Неявная схема Последовательная схема Параллельная схема

1 433 125 66

2 225 67 36

4 142 46 24

8 102 48 25

Таблица 1: Зависимость времени счета (в минутах) от числа процессоров

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

Рис. 5: Разбиение области по процессорам

Основные результаты работы:

1. Предложены и исследованы аддитивные схемы расщепления по направлениям для псевдопараболической модели двойной пористости;

2. Построены безусловно устойчивые схемы расщепления по процессам для общей модели двойной пористости;

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

Список литературы

1. Вабищевнч, П. Н., Григорьев, А. В. Схемы расщепления для псевдопараболических уравнений // Дифференциальные уравнения. 2013. Т. 49, №7. с. 837-843.

2. Григорьев, А. В. Численное моделирование фильтрации в трещиновато-пористой среде // Математические заметки ЯГУ. 2013. Т. 20, выпуск 2. с. 237-245.

3. Григорьев А. В. Preparing input data for modeling the development of oil fields // Труды международных конференций по математическому моделированию, Под редакцией В.И. Васильева. Якутск: Сфера. 2012. с. 31-34.

4. Григорьев А. В. Массивно-параллельное решение уравнения Пуассона с использованием технологии CUDA // Труды международных конференций по математическому моделированию, Под редакцией В.И. Васильева. Якутск: Сфера. 2012. с. 308-314.

5. Gaspar F. J., Grigoriev A. V., Vabishchevich P. N. Explicit-implicit splitting schemes for some systems of evolutionary equations // International Journal of Numerical Analysis & Modeling. Accepted for publication, 2013.

6. Григорьев А. В. Реализация модели двойной пористости с применением гетерогенных вычислений CUDA // Сборник тезисов, Математическое моделирование развития северных территорий 2012, с. 139-140. 2012.

7. Григорьев А. В. Modelling of filtration in fractured porous media // Тезисы Международной конференции: Суперкомпьютерные технологии математического моделирования Якутск: Издательский дом СВФУ, с. 32. 2013.

8. Васильева М. В., Григорьев А. В., Захаров П. Е. Параллельное программирование на основе библиотек: учебное пособие // Якутск: Издательско-полиграфический комплекс СВФУ. 2011, -с. 93.

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ФИЛЬТРАЦИИ В ТРЕЩИНОВАТО-ПОРИСТЫХ СРЕДАХ НА ОСНОВЕ МОДЕЛИ ДВОЙНОЙ ПОРИСТОСТИ

автореферат

Подписано в печать 18.11.13. Формат 60x84/16. Гарнитура «Тайме». Печать офсетная. Печ. л. 1,25. Уч.-изд. л. 1,5. Тираж 100 экз. Заказ X» 358 Издательский дом Северо-Восточного федерального университета, 677891, г. Якутск, ул. Петровского, 5.

Отпечатано в типографии ИД СВФУ

Министерство образования и науки Российской Федерации Северо-Восточный федеральный университет им М. К. Аммосова

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

04201452099

Григорьев Александр Виссарионович

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ФИЛЬТРАЦИИ В ТРЕЩИНОВАТО-ПОРИСТЫХ СРЕДАХ НА ОСНОВЕ МОДЕЛИ ДВОЙНОЙ ПОРИСТОСТИ

Специальность: 05.13.18 — Математическое моделирование, численные методы и комплексы программ

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

Научный руководитель: д.ф.-м.н, профессор Вабищевич Петр Николаевич

Содержание

Введение 5

1 Моделирование двойной пористости в псевдо-параболическом приближении 14

1.1 Схема с весами для псевдо-параболической модели двойной пористости ..................................................................14

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

1.1.2 Разностная схема................................................19

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

1.2 Реализация псевдо-параболической задачи на основе метода конечных элементов............................28

1.2.1 Конечно-элементная реализация................28

1.2.2 Модельная задача........................30

1.2.3 Вычислительный алгоритм...................31

1.2.4 Результаты расчетов.......................32

1.3 Схемы расщепления для псевдо-параболической модели двойной пористости................................48

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

1.3.2 Векторная задача........................52

1.3.3 Аддитивные векторные схемы.................54

1.3.4 Численная реализация.....................59

2 Моделирование двойной пористости на основе системы уравнений 60

2.1 Модель Баренблатта...........................60

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

2.1.2 Постановка обезразмеренной двумерной задачи.......63

2.1.3 Вычислительный алгоритм...................64

2.1.4 Результаты расчетов.......................65

2.2 Схемы расщепления...........................71

2.2.1 Последовательная схема расщепления............71

2.2.2 Параллельная схема расщепления...............72

2.2.3 Более общие задачи.......................73

2.3 Явно-неявные схемы для систем уравнений.............76

2.3.1 Введение.............................76

2.3.2 Начально-краевые задачи для систем уравнений ......79

2.3.3 Схема с весами.........................82

2.3.4 Схемы с диагональным оператором..............84

2.3.5 Общий случай..........................90

3 Моделирование процесса фильтрации в ненасыщенном грунте с применением модели двойной пористости 91

3.1 Модельная задача просачивания воды в грунт............92

3.1.1 Введение.............................92

3.1.2 Уравнение Ричардса ......................93

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

3.1.4 Вычислительный алгоритм...................97

3.1.5 Результаты расчетов.......................98

3.2 Трехмерная задача просачивания воды в грунт............111

3.2.1 Введение.............................111

3.2.2 Результаты расчетов.......................112

3.3 Применение схем расщепления для трехмерной задачи.......121

3.3.1 Введение.............................121

3.3.2 Схемы расщепления ......................122

Заключение 125

Литература 126

Введение

Исследование современных прикладных проблем базируется на применении информационных технологий, интеллектуальным ядром которых выступает математическое моделирование [3,6,29,47,52,54,65,83]. Сформировалась новая технология научных исследований — вычислительный эксперимент. В рамках триады A.A. Самарского модель — алгоритм — программа исследуются важнейшие научно-технические проблемы современности. Вычислительные средства (компьютеры и численные методы) делают возможным описание свойств исследуемого объекта с необходимой полнотой и детальностью на основе адекватных математических моделей, которые включают системы связанных друг с другом нестационарных нелинейных уравнений с частными производными, системы обыкновенных дифференциальных и алгебраических уравнений.

При решении прикладных проблем мы имеем дело с краевыми задачами для систем нестационарных уравнений в частных производных. При построении вычислительных алгоритмов для таких задач решаются проблемы аппроксимации уравнений с учетом соответствующих начальных и граничных условий. Аппроксимация по пространству проводится на основе разностных методов, метода конечных элементов, метода конечных объемов [20,55,60,75]. В настоящее время необходимо решать задачи в сложных расчетных областях. В силу этого метод конечных элементов рассматривается как основная технология для проведения инженерных и научных исследований.

Особые требования предъявляются к выбору аппроксимаций по времени при численном решении задач для систем уравнений [35, 56, 61]. Помимо общих

условий аппроксимации и устойчивости необходимо иметь в виду и вопросы вычислительной реализации построенных схем — решения соответствующей сеточной задачи на новом временном слое. В этом плане наиболее обнадеживающие результаты связываются с построением специальных аддитивных оператор-но-разностных схем (схем расщепления) [17,24].

При исследовании разностных схем для нестационарных задач математической физики широко используется общая теория устойчивости (корректности) операторно-разностных схем [9,20,28,77]. В настоящее время получены точные (совпадающие необходимые и достаточные) условия широкого класса двух- и трехслойных разностных схем в конечномерных гильбертовых пространствах. Необходимо особо подчеркнуть конструктивность общей теории устойчивости операторно-разностных схем, в которой критерии устойчивости формулируются в виде легко проверяемых неравенств для операторов. Среди наиболее важных обобщений отметим использование общей теории устойчивости для некорректных эволюционных задач [21,23,78,80] и для исследования проекционно-разностных схем (схем конечных элементов) [9,22,79]. Получены также новые априорные оценки устойчивости в интегральных по времени нормах [26], на основе которых исследуется, в частности, сходимость разностных схем для задач с обобщенными решениями. Особого внимания заслуживают полученные априорные оценки сильной (коэффициентной) устойчивости при различных предположениях о возмущении операторов (коэффициентов) дифференциальной и разностной задач [25,59].

Аддитивные схемы (схемы расщепления) применяются при решении различных нестационарных задач [17,20,24,34] и предназначены для более эффективной вычислительной реализации построенных схем, для нахождения решения соответствующей сеточной задачи на новом временном слое. Переход к цепочке более простых задач позволяет построить экономичные разностные схемы — расщепление по пространственным переменным (локально-одномерные схемы). В ряде случаев полезно отделить подзадачи различной природы — расщепление

по физическим процессам. Регионально-аддитивные схемы (схемы декомпозиции области) ориентированы на построение вычислительных алгоритмов для параллельных компьютеров.

Основные теоретические результаты об устойчивости и сходимости аддитивных схем получены для скалярных эволюционных уравнений первого порядка и, в некоторых случаях, для уравнений второго порядка [17,20,24,34]. Для вычислительной практики значительный интерес представляют схемы расщепления для систем эволюционных уравнений. Например, для векторных задач отдельные компоненты вектора неизвестных связаны друг с другом. В этом случае использование тех или иных схем расщепления ориентировано на то, чтобы получить хорошие задачи для отдельных компонент решения на новом временном слое.

Для стандартных параболических и гиперболических систем уравнений с самосопряженным эллиптическим оператором аддитивные схемы построены в [20] на основе принципа регуляризации разностных схем. Для систем уравнений схемы расщепления могут строиться с использованием треугольного расщепления оператора задачи на сумму сопряженных друг другу операторов — попеременно-треугольного метода A.A. Самарского. Такие аддитивные схемы используются в работе [62] для динамических задач упругости. Аналогичный подход применяется [11,76] для задач несжимаемой жидкости с переменной вязкостью. Аддитивные схемы для нестационарных задач электродинамики рассмотрены в работе [8].

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

Различные классы аддитивных схем построены для векторных задач [76]. В случае, когда отдельные компоненты вектора неизвестных связаны друг с другом, использование тех или иных схем расщепления ориентировано на то, чтобы получить хорошие задачи для отдельных компонент решения на новом временном слое. Для параболических и гиперболических систем уравнений с самосопряженным эллиптическим оператором локально-одномерные аддитивные схемы построены в [20] на основе принципа регуляризации разностных схем. Для систем уравнений эффективные схемы расщепления могут строиться на основе попеременно-треугольного метода A.A. Самарского, который обычно рассматривается как итерационный метод [20,30]. Такой подход реализован, в частности, в [62] для динамических задач упругости, в [11] — для задач несжимаемой жидкости с переменной вязкостью. Аддитивные схемы для нестационарных векторных уравнений первого и второго порядка, типичных для задач электродинамики, построены в работе [8].

Аддитивные операторно-разностные схемы для систем эволюционных уравнений строятся при зацеплении операторов по пространству. В некоторых случаях имеет место зацепление по производным компонент вектора решения по времени. Поэтому представляет интерес построение аддитивных операторно-разностных схем с расщеплением оператора при производной по времени. Теория и практика построения такого рода схем расщепления в настоящее время находится в зачаточном состоянии. Фактически первой работой для схем расщепления для задач с аддитивным представлением оператора при производной по времени для эволюционных уравнений первого порядка является статья [89]. В ней предложены и исследованы новые векторные аддитивно-операторные схемы при расщеплении оператора при производной по времени на сумму положительно определенных самосопряженных операторов. Такого рода схемы нельзя напрямую использовать для систем эволюционных уравнений при зацеплении по производным по времени.

Большое прикладное значение в строительстве гидротехнических сооружений, в мелиорации, водоснабжении, при добыче нефти и газа имеют исследования движение жидкости (воды, нефти) или газа (воздуха, природного газа) сквозь пористую среду. Математическая модель динамики жидкости в пористой среде включает дифференциальные уравнения, которые в той или иной форме выражают законы сохранения массы и количества движения [41,45]. Прежде всего, используются уравнения неразрывности, которые соответствуют закону сохранения массы для каждой отдельной фазы. В пористой среде уравнения движения записываются в виде закона Дарси, который связывает скорость с давлением.

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

Математическое моделирование течений многофазной жидкости в пористых средах имеет важное прикладное значение при добыче нефти и газа. Традиционно гидродинамические симуляторы строятся на основе трехфазной модели (three-phase black oil) [37, 69]. Эти прикладные математические модели являются существенно нелинейными и трудными для исследования [87,94]. Вторая особенность математических моделей течений многофазных жидкостей, которая характерна и для линейных задач, проявляется в замыкании системы уравнений на основе постоянства суммы всех насыщенностей. Такие алгебраические составляющие модели необходимо учитывать при построении вычислительных алгоритмов решения задач фильтрации многофазной жидкости [43,50].

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

проблемах мы должны учитывать, что в породах, кроме пор, имеется развитая система трещин. Математические модели движения жидкости в такой среде были разработаны в конце 50-х годов Г.И.Баренблаттом, Ю.П.Желтовым, И.Н.Кочиной [4,5]. В современной литературе эта модель двойной пористости (в порах и в трещинах) известна как модель Баренблатта, а в качестве основной цитируемой работы выступает [39] — английский вариант статьи [5]. Модель характеризуется наличием обмена давлениями между фазами.

Подобная математическая модель с перетоком строится для моделирования фильтрации жидкости в многопластовых системах. Моделирование фильтрации жидкости для практических нужд основывается чаще всего на использовании площадных моделей (фильтрация в плане) [14,18]. В этом случае используются приближения с осреднением фильтрационных потоков по толщине водоносного слоя. Типичной является ситуация, когда отдельные водоносные слои перемежаются слабопроницаемыми пластами. Математические модели фильтрации в многопластовых системах строятся на основе предположения (модель Митяева-Гиринского) о преимущественном продольном течении жидкости в водоносных слоях и поперечном течении в разделяющих слоях. Эти математические модели могут быть обоснованы на основе теории гомогенизации [31].

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

приводит к одному уравнению, которое является псевдопараболическим уравнением [81,82]. Хорошо известны схемы расщепления по пространственным переменным для параболических задач. Интересно построить схемы расщепления по направлениям для краевых задач для псевдопараболических уравнений.

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

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

Для инженерных и научных вычислений широко используется библиотека PET Sc (Portable Extensible Toolkit for Scientific Computation). Этот программный инструментарий поддерживает современные парадигмы параллельного программирования на основе стандарта MPI. Основное внимание в PETSc уделяется численному решению линейных и нелинейных систем уравнений, которые возникают при приближенном решении краевых задач для уравнений с частными производными.

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

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

В первой главе ра