Содержание к диссертации
Введение
Глава 1. Модифицированная модель для описания фильтрации однородной жидкости 18
1.1 Постановка задачи 19
1.2 Введение стабилизационной поправки в уравнение неразрывности 22
1.3 Моделирование одномерной задачи фильтрации однородной жидкости 26
1.4 Описание полученных результатов 31
Глава 2. Обобщение модели на двумерный случай. Использование вложенных сеток 33
2.1 Постановка задачи и разностная схема 33
2.2 Результаты численных расчетов 36
2.3 Реализация метода вложенных сеток 37
Глава 3. Численное моделирование трехмерных задач фильтрации на вложенных сетках 41
3.1 Постановка задачи для трехмерного случая 41
3.2 Результаты расчетов 45
Глава 4. Численное моделирование неустойчивости Саффмана - Тейлора в пористой среде. Образование "языков вытеснения" 48
4.1 Линейный анализ 49
4.2 Постановка задачи численного моделирования неустойчивости Саффмана - Тейлора 55
4.3 Результаты численных расчетов 59
Заключение 61
Список литературы
- Введение стабилизационной поправки в уравнение неразрывности
- Результаты численных расчетов
- Постановка задачи для трехмерного случая
- Постановка задачи численного моделирования неустойчивости Саффмана - Тейлора
Введение к работе
Процессы фильтрации широко проявляют себя в самых различных направлениях хозяйственной деятельности человека. Фильтрация связана с движением жидкостей или газов в пористых средах [1-4], которые могут быть как природного происхождения: нефтяные или газовые месторождения, водоносные пласты, питающие реки и озера, так и искусственные очистительные фильтры, химические реакторы. Все шире применяются модели фильтрации и в медицине, когда исследуются течения крови во внутренних органах человека [5].
Последовательное изучение процессов фильтрации началось с подземной гидрогазодинамики - науки о движении жидкостей, газов и их смесей в пористых и трещиноватых средах. По своей сути она является одним из специальных разделов общего курса механики жидкостей. С другой стороны, подземная гидрогазодинамика является теоретической базой для описания процессов фильтрации при разработке нефтяных и газовых месторождений и обеспечивает решение широкого круга прикладных задач в практической деятельности специалистов-нефтяников. Подземная гидрогазодинамика при решении стоящих перед ней теоретических и практических задач пользуется всеми известными в механике жидкостей приемами и методами. Объектом изучения подземной гидрогазодинамики является поток жидкости и газа в пористой среде, называемый фильтрационным потоком; движущиеся в пласте жидкости и газы рассматриваются как непрерывные сплошные среды, на которые распространяются свойства, присущие сплошным средам, и законы механики сплошных сред.
Теория фильтрации имеет глубокие исторические корни (см. очерк развития теории фильтрации в [6]). В 1852 - 1855 гг., производя опыты по фильтрации в песчаных грунтах, французский инженер Г. Дарси установил линейную связь между скоростью фильтрации воды и потерями напора, называемую законом фильтрации, или законом Дарси [7]. Т.е. во второй половине XIX столетия началось развитие теории фильтрации.
Первые теоретические исследования фильтрации, основанные на линейном законе фильтрации, были начаты, Ж. Дюпюи [8]. Они касались случаев фильтрации в грунтовых руслах и осесимметричной фильтрации к колодцам, доходящим до горизонтальной водоупорной поверхности. Эти исследования относятся к 60-м годам XIX столетия.
Позднее, начиная с 80-х годов XIX века, Ф. Форхгеймер рассмотрел более сложные задачи по фильтрации при наличии горизонтального напора. Однако общей теории и общих дифференциальных уравнений фильтрации до 1889 г. не было. Первая работа в этом направлении была написана выдающимся русским математиком и механиком Н. Е. Жуковским. Работа эта называлась "Теоретическое исследование о движении подпочвенных вод". В ней Жуковский вывел дифференциальные уравнения фильтрации [9].
В 1922 г. теория фильтрации получила новый толчок в своём развитии благодаря работе Н. Н. Павловского "Теория движения грунтовых вод под гидротехническими сооружениями и её основные приложения". Этот труд послужил фундаментом, на котором развивалось гидротехническое направление советской школы фильтрации.
Также, в 20-х годах XX века в результате теоретических и экспериментальных исследований появился ряд работ Л. С. Лейбензона, посвященных фильтрации нефти, газа и газированной жидкости. В 1934 г. была опубликована его монография, посвященная теории фильтрации в связи с проблемами добычи нефти и газа. Именно ему принадлежит приоритет в постановке и решении ряда задач нефтегазовой и подземной гидромеханики [10,11].
Выдающийся вклад в теорию фильтрации в нефтегазоводоносных пластах внесли академик С. А. Христианович [12,13], профессоры Б. Б. Лапук [14], И. А. Чарный [15,16], В. Н. Щелкачёв [17 - 20]. Написанные ими монографии и учебники стали классическими, основополагающими. Они имеют большое научно-методологическое значение.
В послевоенный период теория фильтрации развивается трудами учёных, среди которых следует отметить работы М. Т. Абасова [21, 22], М. Г. Алишаева [23], И. М. Аметова, Е. Ф. Афанасьева, Г. И. Баранблатта [24], 10. П. Борисова [25, 26], С. Н. Бузинова [27], В.Я.Булыгина [28], Г. Г. Вахитова [29], М. М. Глоговского, Г. Л. Говоровой, А. Т. Горбунова [30], М. А. Гусейн-Заде [31, 32], В.Л.Данилова [33], Ю.В.Желтова [34], Ю.П.Желтова [35, 36], С.Н. Закирова [37], Г.А.Зотова [38], В. М. Ентова [4], Р.Г.Исаева, Ю. П. Коротаева [39], А. К. Курбанова, Е. М. Минского, Ю. М. Молоковича [40], А. X. Мирзаджанзаде [41,42], Н. Н. Непримерова [43], В. Н. Николаевского [44, 45], Г. Б. Пыхачева [46], Г.В.Рассохина, М. Д. Розенберга [47], Е.С.Ромма [48], Э. В. Соколовского [49], М. Л. Сургучева [50], М. М. Саттарова, Ф. А. Требина [51], Э. Б. Чекалюка [52], М. В. Филинова, М. И. Швидлера [53], И. Д. Умрихина, М. М. Хасанова [54], А. Л. Хейна, Д. А. Эфроса [55] и др.
Работы этих учёных и их учеников обеспечили успешное развитие подземной гидромеханики - теоретической основы теории и практики разработки нефтяных и газовых месторождений, что способствовало ускоренному развитию нефтегазодобывающей промышленности нашей страны.
Широкие исследования в области нефтегазовой подземной гидромеханики проводятся также и за рубежом. Большое значение, для развития технологии нефтеотдачи за рубежом, имеют работы, по теории фильтрации однородных и неоднородных жидкостей, крупнейшего американского специалиста М. Маскета, который хорошо известен нам, по своим двум капитальным монографиям [56, 57]. Большое значение для развития теории фильтрации в целом имеют фундаментальные работы А. Э. Шейдеггера (Канада) [58], Н. Кристеа (Румыния) [59], Р. Коллинза (США) [60], X. Азиза и Э. Сеттари (США) [61, 62]. В гидротехническом аспекте задачи двухфазной фильтрации впервые были рассмотрены П. Я. Полубариновой-Кочиной [63] при изучении плоской напорной неустановившейся фильтрации пресной и соленой воды в основании бетонной плотины гидроузла. Без учета силы тяжести двухфазная фильтрация для случая прямолинейно-параллельного вытеснения рассматривалась С. Баклеем и М. Левереттом, а также существенное значение имеют основы теории многофазной фильтрации, также разработанные С. Баклеем и М. Левереттом [64,65], а позже независимо от них А. М. Пирвердяном [66].
Благодаря полученным результатам в области фундаментальной науки, стало возможным плодотворное применение гидродинамических расчетов в гидрогеологических изысканиях, прогнозах подпора грунтовых вод в районах гидросооружений и орошаемых территорий; оценках эксплуатационных запасов подземных вод и т.п.
С появлением в 60-е годы мощных ЭВМ в решении практических задач фильтрации все большее место занимает численное моделирование. Однако аналитические методы и в этом случае остаются востребованными не только для разработки и тестирования численных алгоритмов, но и для глубокого понимания физики явлений, параметрического анализа сложных схем, оптимизации и оценок характеристик фильтрационных полей, в том числе в ситуациях, характеризующихся высокой степенью неопределенности относительно параметров пористой среды, механизмов взаимодействия флюида и матрицы, граничных условий и самих границ области течения. В обзоре [67] охватываются работы по кругу проблем, непосредственно связанные с исследованиями по динамике подземных вод в областях, границы которых известны не полностью. В математической постановке эти проблемы сводятся к краевым задачам для дифференциальных уравнений эллиптического типа в областях с неизвестными границами, которые отыскиваются по заданным на них краевым условиям.
Как показывает практика последних десятилетий, математическое моделирование существенно влияет на развитие большинства направлений современной науки и техники. Задачи, связанные с фильтрацией, в этом смысле, - не исключение. Здесь нужно отметить несколько аспектов. Во-первых, необходимо совершенствование математических моделей теории фильтрации. Во-вторых, необходимо развитие и совершенствование методов реализации этих моделей. Наконец, необходимо учитывать быстрое развитие высокопроизводительных параллельных вычислительных машин, что позволяет решать задачи повышенной сложности, но требует специальных методов и алгоритмов. В 60-х годах выдающимся ученым Н. Н. Яненко [68] было предложено в вычислительной математике создание метода дробных шагов. Метод дробных шагов - это метод построения экономичных (в смысле числа операций) конечно-разностных схем для решения дифференциальных уравнений. Решение задач гидродинамики в сложных областях потребовало развития новых методов построения расчетных сеток, автоматически адаптирующихся к потоку, что позволило, в конечном счете, на порядок повысить точность расчетов.
Математическая модель - это система дифференциальных уравнений, описывающая процесс фильтрации в рассматриваемом конкретном объекте разработки, с заданными начальными и граничными условиями, обеспечивающими единственность решения поставленной задачи.
Этапами подготовки математической модели в данном случае являются:
- создание геологической модели;
- обоснование размерности модели и выбор основных уравнений для описания процесса;
- задание начальных и граничных условий.
В модельных задачах подземной гидрогазодинамики обычно рассматривают упрощенные геологические модели, в которых участок разработки или залежь в целом схематизируется в прямоугольную или круговую область фильтрации с постоянной толщиной пласта. В геологической модели должны быть также заданы основные геолого-физические параметры пластовой системы (пористой среды и фильтрующихся флюидов). В упрощенных моделях - это средние значения этих параметров в моделируемой области.
Для простейших линейно-параллельного и радиального потоков в пласте постоянной толщины применимы одномерные модели фильтрации. Для сложных течений в областях, содержащих произвольные системы скважин, используются двумерные (в случае тонких пластов) и трехмерные (в наиболее общем случае) модели фильтрации.
К основным уравнениям математической модели относятся:
- уравнения неразрывности (законы сохранения массы) для каждой фильтрующейся фазы;
- уравнения движения (обобщенный закон Дарси) для каждой фазы;
-уравнение сохранения энергии (в случае неизотермической фильтрации);
- уравнения состояния;
- дополнительные соотношения, устанавливающие взаимосвязи между фазовыми насыщенностями, между фазовыми и капиллярными давлениями.
Одной из важнейших индустриальных задач является проблема фильтрации нефти. Пример решения таких задач можно найти в работе [69], а также в книге [70].
Основы современных подходов математического моделирования изложены в работах [71 - 89].
Среди множества работ по математическому моделированию процессов фильтрации, опубликованных в последнее десятилетие, хотелось бы отметить работы [90-101].
В последнее время бурно развивается многопроцессорная вычислительная техника. Многие алгоритмы решения задач описываемых уравнениями математической физики, приводят к большим затратам вычислительных мощностей (численное моделирование задач фильтрации относится к этому кругу), и требуют адаптации к архитектуре новой многопроцессорной вычислительной техники. Некоторые общие требования к алгоритмам и задачам для многопроцессорных вычислительных систем отмечены в монографии [102].
Одним из возможных путей развития алгоритмов для многопроцессорной вычислительной техники, как было указано в [102], является разработка явных схем [103] с использованием специальных «регуляризаторов» [104]. При этом в алгоритм или непосредственно в математическую модель вводятся дополнительные члены, которые демпфируют счетные осцилляции. Данные «регуляризаторы» не должны превращать изучаемую задачу в среду, в которой гасятся все естественные возмущения. С другой стороны они должны успешно исключать неустойчивости именно счетного характера. Примером построения таких алгоритмов могут служить кинетически-согласованные разностные схемы, используемые для моделирования течений вязкого газа [102]. Заметим, что для создания подобного рода методов необходимо опираться не только на знание теории разностных схем, но и активно использовать разнообразную физическую информацию о свойствах изучаемого явления.
Важным свойством алгоритма, позволяющим добиться успешной адаптации на многопроцессорные системы [78], является его логическая простота. Это требование наиболее важно в связи с продолжающимся прогрессом многопроцессорной вычислительной техники, сопровождаемым неизбежным в той или иной степени изменением математического обеспечения.
Указанные выше требования к алгоритмам для многопроцессорных вычислительных систем учитывались автором при разработке новых моделей фильтрации с использованием специальных «регуляризаторов» и вложенных сеток.
Целью диссертационной работы является создание алгоритмов решения, численных методов и комплексов программ для решения актуальных задач фильтрации. Численные алгоритмы должны удовлетворять следующим требованиям: 1) использование явных однородных схем, что должно в будущем обеспечить высокую эффективность реализации данных алгоритмов на параллельных вычислительных системах; 2) разрешение больших градиентов течения вблизи скважин, для чего было предложено использовать систему вложенных сеток; 3) возможность решения двумерных и трехмерных задач; 4) переход к многофазной фильтрации.
В качестве первой задачи исследования была поставлена одномерная задача о центрально-симметричной фильтрации. Необходимо было найти численное решение на основе использования простых симметричных явных или полунеявных схем с возможной последующей регуляризацией решения. В качестве второй задачи исследования была поставлена двумерная задача о фильтрации вблизи скважины. Необходимо было описать на сетке скважину и градиенты численного решения, которые велики вблизи скважины. Для разрешения больших градиентов был выбран один из способов «сгущения» сетки вблизи скважины: система вложенных сеток. Следующей задачей была поставлена переформулировка алгоритмов и программ на трехмерный случай с учетом гравитации. Наконец, последняя задача была связана с исследованием устойчивости фильтрационного течения, когда менее вязкая жидкость проталкивает вязкую через пористую среду. При построении математической модели и разностных схем использовался опыт построения кинетически-согласованных разностных схем (КСРС) и квазигазодинамической системы уравнений [102]. При отладке программ проверялась сходимость численных решений при измельчении сетки. Полученные решения сравнивались с точными решениями в тех случаях, когда оно известно [105]. Сравнивались численные решения, полученные разными методами. При разработке алгоритмов и программ фильтрации на системе вложенных сеток учитывался опыт коллег, численно решавших задачу о взрыве сверхновой на системе вложенных сеток [106,107].
В ходе выполнения данной работы создан новый модифицированный алгоритм численного решения задачи фильтрации однородной жидкости, полученный на основе введения регуляризирующей поправки. Данный алгоритм подразумевает использование вложенных сеток, которые в свою очередь служат для разрешения больших градиентов (например, вблизи области скважины). Такая технология вложенных сеток имеет определенные преимущества по сравнению с другими технологиями измельчения, например: 1) простота расчетных схем из-за регулярности сеток, 2) существенная экономия ресурсов - четырем шагам на мелкой сетке соответствует только один шаг на более крупной сетке. Выполнена реализация разработанных методов и для трехмерных задач с учетом гравитации.
Разработанные методы, алгоритмы и комплексы программ могут быть использованы, как для фундаментальных исследований задач однофазной фильтрации, так и в приложениях, например, для расчетов водоносных пластов задач фильтрации нефти. Разработанный комплекс программ позволяет рассчитывать трехмерную область фильтрации с одной скважиной, разрешаемой системой вложенных сеток. Комплекс может быть обобщен на произвольное количество скважин. Разработанный метод с регуляризирующей поправкой может быть использован для больших скоростей фильтрующего вещества (закон Форхаймера вместо закона Дарси). Реализованная методика вложенных сеток дает простой и эффективный механизм расчета задач гидромеханики и фильтрации, в которых есть выделенная точка (ось) и градиенты сильно зависят от расстояния до точки (оси). Выполнено моделирование неустойчивости Саффмана - Тейлора.
Результаты работы устно докладывались и обсуждались на конференциях и семинарах:
- Японско - Российский семинар "Турбулентность, неустойчивости и параллельные вычисления" (г. Осака, Япония, 22 - 24 октября, 2001);
- XXXXV научная конференция МФТИ "Современные проблемы фундаментальных и прикладных наук" (г. Долгопрудный, МФТИ, 29-30 ноября, 2002);
- VIII Японско - Российский симпозиум "Вычислительная гидрогазодинамика" (г. Сендай, Япония, 24 - 26 сентября, 2003);
- Японско - Российский семинар "Турбулентность и неустойчивости" (г. Токио, Япония, 29 - 30 сентября, 2003);
- Российско - Индийская международная конференция "Высокопроизводительные вычисления в науке и промышленности" (г. Москва, 16-20 июня, 2003);
- Российско - Японская международная конференция "Турбулентность и неустойчивости" (г. Москва, 21-24 сентября, 2004).
Основное содержание диссертации опубликовано в работах автора [108 — 110] и тезисах докладов [111-113].
Диссертация состоит из введения, четырех глав, заключения и списка литературы. Объем работы составляет ПО страниц. Работа содержит 35 рисунков и список цитированной литературы из 166 наименований. В первой главе, состоящей из четырех разделов, описывается новая модель фильтрации однородной жидкости на основе закона Дарси. При выводе модели используется опыт построения кинетически-согласованных разностных схем и квазигазодинамической системы уравнений.
Раздел 1.1 посвящен описанию математической модели фильтрации и постановки задачи, основанной на законе Дарси. Закон Дарси дополняется классическим уравнением неразрывности и уравнением состояния, которое связывает давление и плотность фильтрующейся жидкости. Для полноты постановки задачи задаются граничные условия и начальные данные.
В разделе 1.2 представлено модифицированное уравнение неразрывности в пространственно-многомерном случае, которое создано на базе уравнения неразрывности путем добавления «регуляризатора» к правой части уравнения. Впоследствии, правую часть этого уравнения будем называть «регуляризатором», или стабилизационной поправкой, или диссипативным членом.
Численная методика, постановка одномерной тестовой задачи и результаты тестовых одномерных расчетов на основе новой модели фильтрации приводятся в разделе 1.3. В одномерном осесимметричном случае жидкость притекает к гидродинамически совершенной скважине.
В разделе 1.4 описываются полученные численные результаты для одномерного осесимметричного случая.
Во второй главе, состоящей из трех разделов, рассматривается двумерный случай. В разделе 2.1 излагается постановка задачи. В разделе 2.2 описываются результаты тестовых двумерных расчетов.
Для эффективного разрешения областей с большими градиентами физических величин был реализован метод вложенных сеток, описанный в разделе 2.3. В зоне с большими градиентами используются более подробные вложенные сетки. Главной проблемой расчета на системе вложенных сеток является передача значений физических величин из одной сетки в другую.
Для апробации реализованного метода вложенных сеток была рассмотрена задача фильтрационного потока вблизи горизонтальной скважины. Течение вблизи скважины должно разрешаться на самой подробной вложенной сетке. Расчетная область покрывается набором вложенных друг в друга сеток с одинаковым числом ячеек, но последовательно уменьшающимся в 2 раза абсолютным размером. Таким образом, расчетная область вблизи скважины имеет гораздо более детальное разрешение по сравнению с периферией.
Для выполнения условия устойчивости необходимо, чтобы четырем временным шагам вычислений на сетке т-го уровня соответствовал один шаг вычислений по времени на сетке (т - l)-ro уровня (если T h2). Граничные значения, необходимые для вычисления величин на внешней границе т-го уровня с ячейками сетки (m +1)-го уровня, получались в результате интерполяции, а граничные значения на внутренней сетке - путем осреднения по ячейкам сетки (ш -1)-го уровня. В вычислениях использовалось до семи уровней вложенных равномерных сеток. Вычисления на системе вложенных сеток реализованы в виде рекурсивной процедуры.
Вложенные сетки, позволяют резко уменьшить компьютерное время вычислений при заданной точности решения. При расчетах использовались различные типы интерполяции для передачи значений физических величин между соседними вложенными сетками. Был проведен сравнительный анализ интерполяций и было показано, что нет необходимости использовать интерполяцию более высокой степени точности, чем линейная. Разработанная методика моделирования фильтрации является явной, что говорит о перспективности ее эффективной реализации на многопроцессорных вычислительных системах.
В третьей главе, состоящей из двух разделов, рассматривается пространственный случай. В разделе 3.1 излагается постановка задачи. Записываются расчетные схемы для трехмерного случая. В разделе 3.2 представлены результаты расчетов вертикальной гидродинамически несовершенной скважины.
В случае, если менее вязкая жидкость вытесняет более вязкую в пористой среде, фронт разделяющий их неустойчив. Линейный анализ и первые эксперименты в ячейке Хеле-Шоу (Hele-Shaw) связывают с именами Саффмана и Тейлора (Saffman - Taylor) [114]. Важнейшей причиной для изучения этой проблемы является то, что она тесно связана с многими технологическими процессами, где присутствует течение в пористой среде. Например, феномен "языков вытеснения" или "пальцев" часто встречается, когда нефть извлекается из под слоя воды. Четвертая глава, состоящая из трех разделов, посвящена моделированию этой неустойчивости. Раздел 4.1 является вводным, в нем аналитически рассматривается линейная стадия развития неустойчивости [115]. В разделе 4.2 рассматривается математическая модель процесса вытеснения жидкостей [33, 100, 101, 116]. Используемая в настоящей работе модель фильтрации включает в себя следующие предположения:
? жидкость считается несжимаемой, пористость однородной;
? выполняется закон Дарси, причем, подвижность зависит от насыщенности;
? насыщенность удовлетворяет конвективно-дисперсионному уравнению. Дисперсия предполагается однородной и изотропной.
С помощью этой модели возможно численное исследование данной неустойчивости и, с учетом линейной зависимости коэффициента подвижности от насыщенности модель сводится к той же системе уравнений и дополнительному уравнению переноса для коэффициента подвижности.
Коэффициент подвижности теперь характеризует соответствующую жидкость и зависит от положения точки в пространстве. В качестве граничных условий задаются давление на внутренней сфере (окружности) вытеснения и внешнее давление (на внешней границе области). Положение фронта вытеснения не выделяется, а вычисляется сквозным образом. Далее в тексте диссертации приводятся постановки задач.
Результаты численного моделирования по вышеописанной модели описываются в разделе 4.3.
В заключении сформулированы основные результаты диссертационной работы:
Разработан оригинальный метод решения задач фильтрации однородной жидкости, основанный на введении диссипативного члена в уравнение неразрывности.
Разработана методика вложенных сеток, которая эффективно позволяет разрешать большие градиенты физических величин вблизи выделенных точек (скважин).
Создан комплекс программ для численного решения двумерных и трехмерных задач фильтрации на вложенных сетках. Выполнены многочисленные тестовые расчеты. Использование регулярных сеток и явных схем предполагает быстрый перенос на многопроцессорные вычислительные системы.
Поставлена и численно рассмотрена задача развития неустойчивости фронта вытеснения жидкости менее вязкой жидкостью.
Автор выражает глубокую благодарность своему научному руководителю за постоянное внимание, поддержку и совместную научную работу: чл.-корр. РАН Четверушкину Борису Николаевичу, д.ф.-м.н. Чечеткину Валерию Михайловичу за ценные обсуждения и полезные советы, а также всему коллективу ИАП РАН, особенно д.ф.-м.н. Трошкину Олегу Валентиновичу и к.ф.-м.н. Семенову Илье Витальевичу, за тесное научное сотрудничество. Отдельно хотелось бы поблагодарить к.ф.-м.н. Опарина Алексея Михайловича за неоценимую помощь в научной работе.
Введение стабилизационной поправки в уравнение неразрывности
Бахвалов Н. С. и Панасенко Г. П. [137] показали, что уравнение Дарси (6а) вытекает из уравнения движения системы уравнений Навье - Стокса в пренебрежении конвективными членами и его осреднении на характерном масштабе /, когда фильтрующую породу можно считать однородной.
Заметим, что рассмотренная выше постановка задачи противоречива по своему построению. В самом деле, если уравнение (6а) описывает процессы с минимальным характерным масштабом /, то уравнение (66) позволяет описывать сложное микроскопическое течение флюида на расстоянии размера зерен породы и меньшем.
Чтобы преодолеть это противоречие, попробуем видоизменить уравнение (66) с учетом характерного масштаба осреднения /. Для этого воспользуемся уже ранее полученным опытом построения кинетически-согласованных разностных схем (КСРС) и квазигазодинамической системой уравнений (КГСУ)[102].
При построении КСРС в явном виде используется тот факт, что одночастинная функция распределения, моментами которой и являются газодинамические параметры, слабо меняется на расстоянии порядка длины свободного пробега молекулы li . Проведем аналогичное построение для уравнения неразрывности, заменяя /, на размер /, на котором уже становится несущественной микроструктура фильтрующей породы.
Рассмотрим функцию распределения f(x, v) для фильтрующейся жидкости, где х - пространственная координата, a v - мгновенная скорость молекулы жидкости. Для моментов функции / естественно потребовать выполнения следующих свойств: J/(x,v)dV = /?(x), J/(X,V)VAV = /7(X)U (X), где р{\) - плотность флюида, а и (х) - мгновенная скорость его макроскопического движения (см. Рис. 1), этот рисунок демонстрирует различие между мгновенной скоростью жидкости, при обтекании зерен и (х) и скоростью фильтрации (поток массы) и. Для описания изменения функции распределения рассмотрим кинетическое уравнение Энского [125]: +, = Ф(х,у). (7)
Его отличие от кинетического уравнения Больцмана состоит в наличии более сложного интеграла столкновения Ф(х,у), в котором для случая газов высокой плотности и жидкостей учтен тот факт, что средний объем пространства, приходящийся на одну молекулу, сопоставим с объемом самой молекулы.
Но, для целей настоящей работы понадобится только закон сохранения массы. Необходимый для его построения момент от интеграла столкновения уравнения Энскога равен нулю: J0(x,v /v = 0, (8) также как и для уравнения Больцмана .
Явная разностная аппроксимация уравнения (7), которую для простоты выпишем в пространственно-одномерном случае, примет вид: 2Ах 2 Ах2 Проинтегрируем уравнение (9) с весом т (где т - масса молекулы) по скоростям v, используя свойство (8) интеграла столкновений уравнения Энскога. При этом, так же как и при выводе квазигазодинамической системы ([125, глП, 5]) заменим v на характерную скорость, в качестве которой удобно взять скорость звука с. При этом получим: (10) А РГ-Р!(У) , (puYM-(pu)lgf/(v)+ti(v)jAxc } J At 2Ах
Заметим, что разностная схема (9) получена при условии кусочно-постоянного по пространству вида функции распределения /: A v) = /(v),xe[x(.,x,.+1]. Заметим, что другие моменты Ф с сумматорными инвариантами v2 Ф(Х V)—UV Ф U уже не равны нулю. Этот факт использовался в работе [138] при построении кинетических J 2 схем для плотных газов и жидкостей.
Зададимся вопросом, связано ли ограничение на Дх снизу только с возможностями вычислительной техники или имеет более глубокий физический смысл? Значение скорости фильтрации и получено путем осреднения по пространству на характерный размер / 102 зерен породы скорости действительного гидродинамического течения и . Поэтому в качестве ограничения снизу на Ах следует взять величину /. Учитывая этот факт [ПО], запишем дифференциальное приближение уравнения (10) в пространственно-многомерном случае, путем добавления «регуляризатора» к правой части уравнения, в виде: Я 7 — + div(/?u) = div—gradp, (11) dt 2 где / 10 -102 зерен породы - обсуждаемый выше характерный масштаб осреднения по пространству, с - характерная скорость распространения возмущения, в качестве которой можно выбрать скорость звука в фильтрующейся жидкости. В последствие, правую часть уравнения (11) будем называть «регуляризатором», стабилизационной поправкой или диссипативным членом.
Уравнение (11) с законом Дарси (1) будем использовать в качестве модели для описания фильтрации жидкости. Обсудим некоторые физические предположения, положенные в основу вывода уравнения (11).
Результаты численных расчетов
Для начала приведём классификацию типов сеток, используемых при решении широкого круга задач: Сетки (Grids) 1. Статические (Static): a. структурированные (structured), b. неструктурированные (unstructured). 2. Адаптивные (Adaptive): a. вложенные (nested), b. движущиеся (moving), c. неструктурированные (unstructured). Эта классификация взята из источников [139 - 142].
В настоящее время проводятся интенсивные исследования по разработке методов построения сеток и по внедрению этих методов в численные алгоритмы решения фундаментальных и прикладных проблем [143 - 146]. Более ранние исследования в этой области также связаны с работой Яненко Н. Н. [68]. Описание некоторых алгебраических, дифференциальных и вариационных методов построения адаптивных сеток дано в монографии [147], обзорах по конструированию разностных [148] и адаптивных [149, 150] сеток, в том числе для сложных геометрий [151], неструктурированных сеток [152].
В книге [153] изложены кратко некоторые современные методы построения адаптивных сеток. И представлены три основные класса сеток, получивших широкое распространение при решении задач в многомерных областях: 1) структурированные, 2) неструктурированные, 3) гибридные.
Для эффективного разрешения областей с большими градиентами физических величин нами был выбран и реализован метод вложенных сеток, который является составной частью метода декомпозиции областей [154]. Забегая вперед, отметим, что данный метод оказался вполне эффективным. В зоне с большими градиентами используются более подробные вложенные сетки. Главной проблемой расчета на системе вложенных сеток является передача значений физических величин из одной области в другую вблизи границ вложенных сеток. Надо сказать, существуют другие способы решения проблемы разномасштабности структур, например, методика прискважинного макроблока [155] и др. [146, 147 - 150,156].
Расчётная область покрывается набором вложенных друг в друга сеток с одинаковым числом ячеек, но последовательно уменьшающимся в 2 раза абсолютным размером ячеек. Сама скважина должна умещаться на самой подробной вложенной сетке (см. Рис. 12), пример с тремя вложенными сетками. На Рис. 13 показан пример с пятью вложенными сетками. Таким образом, расчетная область вблизи скважины имеет детальное разрешение, что приводит к более точному решению.
Для апробации реализованного метода вложенных сеток была рассмотрена задача фильтрационного течения вблизи горизонтальной скважины. В этой задаче существует почти радиальное фильтрационное течение, в котором жидкость притекает к гидродинамически совершенной добывающей скважине. В данном случае, когда мы используем метод вложенных сеток, форма добывающей скважины может быть произвольна (см. Рис. 12).
Для выполнения условия устойчивости необходимо, чтобы четырем временным шагам вычислений на сетке (т)-го уровня соответствовал один шаг вычислений по времени на сетке (т+1)-го уровня (если т-h2). Как раз, Рис. 14 демонстрирует алгоритм вычислений одного шага по времени для трех вложенных сеток. Используемая функция solve(m) проводит вычисления на сетке (т)-го уровня, например, функция solve(O) - проводит вычисления на самой внутренней подробной сетке. Рассмотрим подробнее границу между вложенными сетками (т)-го уровня и (т+1)-го уровня, (см. Рис. 15). Для того чтобы вычислить значение в граничной точке сетки (т)-го уровня (левый черный кружок), необходимо интерполировать значение в дополнительной точке (пустой кружок) по значениям окружающих точек сетки (т+1)-го уровня. Для четных вертикальных индексов, дополнительная точка попадает в центр квадрата, по углам которого находятся узлы сетки (т+1)-го уровня, а для нечетных -дополнительная точка попадает в центр ребра, по краям которого находятся узлы сетки (т+1)-го уровня.
Далее, для того чтобы вычислить значение в граничной точке сетки (т+1)-го уровня (правый верхний черный квадрат), необходимо иметь значение в дополнительной точке (пустой квадрат). Эта точка совпадает с одной из точек сетки (т)-го уровня. Поэтому мы переносим значение с соответствующей внутренней точки сетки (т)-го уровня. В расчетах использовалось до семи уровней вложенных равномерных сеток. Вычисления на системе вложенных сеток реализованы в виде рекурсивной процедуры.
Постановка задачи для трехмерного случая
Кроме так называемых гидродинамически совершенных скважин (т.е. вскрывших пласт на всю толщину и сообщающихся с пластом по всей площади сечения забоя скважины) в реальных ситуациях чаще приходится иметь дело с несовершенными скважинами. К категории несовершенных скважин относятся скважины, вскрывшие пласт не полностью (скважины, несовершенные по степени вскрытия), а также скважины, сообщающиеся с пластом через перфорационные каналы искусственного фильтра (см. [157]). В этом случае моделирование фильтрации необходимо проводить в трехмерной по пространству постановке.
Рассмотрим в качестве трехмерной постановки задачи случай гидродинамически несовершенной по степени вскрытия скважины. Задача рассматривается в декартовых координатах [x,y,z). Ось z направлена вертикально. Область интегрирования представляет собой прямоугольный параллелепипед /) = [-Х,Х]х[-Г,7]х[-7,0]. Плоскости z = -Z и z = 0 ограничивают несущий пласт и непроницаемы для фильтрующейся жидкости. Скважина S с координатами центра (х5 =0,ys = 0) входит в пласт сверху вертикально и доходит до середины пласта. Задан радиус скважины гс = 50. Начальное давление в области интегрирования соответствует гидростатическому P = P0-\g\pz. На верхней границе области оно равно Р0 = 107. В начальный момент включается скважина и давление в ней равно Рс = 106. Скорость в начальный момент нулевая. Поддерживаются давления, равные начальному, на боковых границах области интегрирования и на скважине. Через верхнюю и нижнюю границу пласта поток фильтрационной жидкости отсутствует. Сила тяжести направлена вниз и по абсолютной величине равна g = 103.
Необходимо найти параметры течения внутри области интегрирования по времени вплоть до установления. Для этого введём равномерную по осям х, у и z сетку:
Значения переменных p, P, и, v,w заданы в узлах. Например, обозначение Рщі относится к значению плотности в узле сетки с координатами \xm,yq,zA в момент времени f .
Итак, в трехмерном случае уравнения (1), (11) принимают следующий вид: (39) Как и в двумерной задаче, будем решать уравнения непосредственно в том виде и порядке, в каком они записаны, то есть из закона Дарси (38) найдём компоненты скорости, а из модифицированного уравнения неразрывности (39) - плотность. Давление определяется по плотности из уравнения состояния.
Итак, если известны все компоненты рпЩІІ Рпщ1, ипщП У"Щ1, w"mql на предыдущем временном слое п, то можно вычислить все значения на новом временном слое (w + l) в момент времени /n+1 =t" + т:
Далее, из модифицированного уравнения неразрывности (39) с учетом новых значений скорости можем найти плотность на новом временном слое или давление. Вычисление давления более предпочтительно, поскольку приводит к более точным результатам.
Итак, имеем формулу для вычисления давления из модифицированного уравнения неразрывности (39) с учетом новых значений скорости на новом временном слое:
Постановка задачи численного моделирования неустойчивости Саффмана - Тейлора
Рассматривается простая математическая модель вытеснения вязкой жидкости менее вязкой жидкостью, представляющая собой обобщение модели течения Хеле - Шоу {Hele - Shaw) на трехмерный случай. Уравнения модели для несжимаемой жидкости имеют вид: V(KVP) = 0 U = -KVP a,A:+v( u)=o (65) K = K(x,t), P = P{x,t); K,PER; ueR\ d = 2, 3 Задача рассматривается в области Q с вырезом Qin: Qjn с Q с Rd. Граничные условия:
1. Условие Дирихле для давления на внешней и внутренней границах расчетной области = Р Р =р (66) где Рх Р0 (условие вытеснения). Данное граничное условие соответствует открытой границе в течении Хеле - Шоу.
2. Для уравнения переноса df!C+ V(KU) = 0 на внешней границе ставится условие: к = к0,\1 П 0 дк/дп=0,и-п (67) гп -о
В задаче вытеснения это условие соответствует выполнению к момента, как фронт вытеснения достигает внешней границы.
Начальные данные: = к, до того (68) где Га - положение фронта вытеснения при t = 0. Физическая интерпретация "обобщенной проницаемости" к применительно к задаче о вязком вытеснении в однородной пористой среде и в течении Хеле - Шоу (Hele - Shaw).
Двумерный анализ линейной устойчивости плоского фронта вытеснения, выполненный Саффманом и Тейлором [114], показал развитие неустойчивости при ju01 цх 1. В эксперименте развития неустойчивости Саффмана - Тейлора в течении Хеле - Шоу приводит к образованию "пальцев вытеснения" (целый ряд фотографий структур, образующихся при неустойчивом вытеснении в ячейке Хеле - Шоу, приведен на Рис. 25, 26 (эти рисунки взяты из работы [115])). При распространении фронта вытеснения в пористой среде наблюдается аналогичное явление. В настоящей работе развитие неустойчивости Саффмана - Тейлора в случае двух и трех измерений в течении с гладким начальным контуром рассчитывается на основе модели (65) без явного выделения фронта вытеснения.
Численный метод основан на применении интегральной формы уравнения для "/с" —(\\\кс1а) + \\Kudn = О, (69)
к ячейкам, на которые разбита расчетная область. Здесь da - элемент объема, dn - элемент поверхности расчетной ячейки.
При вычислении интеграла по поверхности ячейки применяется метод Гаусса. Реконструкция решения в точках интегрирования Гаусса на гранях ячеек осуществляется с использованием метода WENO [166], который обеспечивает хорошее разрешение разрывов при численном решении гиперболических уравнений методом сквозного счета. Численные аппроксимации потока на гранях ячеек вычисляются по реконструированным в точках Гаусса "левым" и "правым" значениям с учетом направления скорости. Интегрирование по времени осуществляется с использованием варианта метода Рунге - Кутта: 5 = +А/Ь,(ЛП, (70) C=7C+7 M (V(1)), (71) з ик з " 3 С = r S} +т4? + -A M/?(V(2)). (72) Здесь Ку, интеграл Щ/cdco, и Ьок(ки),Р(1)) аппроксимация интеграла (-jj/aififnj применительно к ячейке с индексами ijk, на временном слое (дробном шаге метода Рунге - Кутта) /. Сеточные функции - к(1), Р{1), представляют решение на слое /, причем, Р{1) определяется на каждом дробном шаге метода Рунге - Кутта как численное решение, соответствующей краевой задачи.
Фронт вытеснения неустойчив при к0/кх \, (см. Рис.27) и устойчив в противном случае Рис. 28. При неустойчивом вытеснении образуется область перемешивания, которая характеризуется взаимопроникновением вытесняющей и вытесняемой жидкостей. Разность скоростей внутренней и внешней границ области перемешивания зависит от отношения подвижностей к0/лг, (или вязкостен), обращаясь в ноль при к0 = к,.
На Рис. 29 и Рис. 30 показаны картины двумерного радиального вытеснения для двух значений к01кх =0.2 и кй1кх =0.447, соответственно. Эффективность вытеснения (отношение объема закачанной жидкости к объему ограниченному внешним фронтом вытеснения) определяется двумя факторами: 1) скоростью расширения области перемешивания, 2) относительным содержанием вытесняющей жидкости в области перемешивания. Первый фактор имеет сильную зависимость от отношения подвижностей (вязкостей), что иллюстрируют Рис. 29, 30. Относительное содержание вытесняющей жидкости по результатам расчетов варьируется в зависимости от размерности задачи. В двумерных расчетах (рассматривалось одномерное и радиальное вытеснение) оно составляло приблизительно 0.5.
