Содержание к диссертации
Введение
Глава 1. Численное моделирование задач пороупругости . 10
1.1. Введение 10
1.2. Постановка задачи 12
1.3. Вычислительные алгоритмы решения задачи пороупругости . 15
1.4. Схемы расщепления для задачи пороупругости 25
1.5. Численные эксперименты 34
Глава 2. Численное моделирование задачи пороупругости для пластин 58
2.1. Введение 58
2.2. Математическая модель 60
2.3. Численное решение задачи прогиба пластин 67
2.4. Задача пороупругости пластин 72
Глава 3. Задачи пороупругости при разработке месторождений нефти и газа 87
3.1. Введение 87
3.2. Математическая модель 89
3.3. Вычислительный алгоритм 93
3.4. Численные эксперименты 94
Заключение 109
Список литературы
- Вычислительные алгоритмы решения задачи пороупругости
- Численные эксперименты
- Численное решение задачи прогиба пластин
- Вычислительный алгоритм
Вычислительные алгоритмы решения задачи пороупругости
Отмеченные выше различные классы аддитивных операторно-разност-ных схем для эволюционных уравнений строятся при аддитивном расщеплении основного оператора (связанного с решением) на несколько слагаемых. В некоторых прикладных задачах интерес представляют задачи при расщеплении оператора при производной по времени. В [140] предложены и исследованы векторные аддитивно-операторные схемы при расщеплении оператора при производной по времени на сумму положительно определенных самосопряженных операторов. В работе [6] выделен класс задач с расщеплением как основного оператора задачи, так и оператора при производных по времени. Для задач, которые возникают при исследовании краевых задач для псевдопараболических уравнений предложены [7] безусловно устойчивые векторные схемы расщепления по направлениям. Явно-неявные схемы для систем эволюционных уравнений рассмотрены в работе [72].
При построении схем расщепления для нестационарных задач пороупру-гости мы ориентируемся на расщепление по физическим процессам. Переход на новый временной слой в этом случае базируется на решении отдельных задач для смещений пористой среды и для давления в порах. В настоящее время некоторые схемы расщепления известны (смотри, например, работы [36, 89–91, 106].), однако их построение базируется на эвристических соображениях инженерного плана, устойчивость таких схем исследуется экспериментально, строгих математических результатов об их устойчивости практически нет. Поэтому актуальной проблемой вычислительной математики является построение абсолютно устойчивых схем расщепления по физическим процессам для задач пороупругости (фильтрационной консолидации), разработка на их основе вычислительных алгоритмов приближенного решения нестационарных задач пороупругости и их применение для решения прикладных проблем фильтрации. Отдельного внимания заслуживают задачи фильтрации для пористых пластин.
Проблемам построения, обоснования и использования схем расщепления по физическим процессам (перемещение и давление) посвящено настоящее исследование. В первой главе обсуждаются проблемы численного решения краевых задач пороупругости (фильтрационной консолидации). Аналогичные математические модели используются при исследовании термоупругости. Базовая система уравнений включает стационарные уравнения Ламе для перемещений и нестационарные уравнения для давления или температуры. Вычислительный алгоритм основан на конечно-элементной аппроксимации по пространству. Для двухслойных схем с весами формулируются стандартные условия устойчивости. Вычислительная реализация таких схем основана на решении системы связанных уравнений для перемещений и давления (температуры). Строятся схемы расщепления по физическим процессам, когда переход на новый временной слой связывается с решением отдельных эллиптических задач для искомых перемещений и давления (температуры). Безусловно устойчивые аддитивные схемы строятся на основе выбора веса трехслойной схемы. Работоспособность предложенных схем расщепления продемонстрирована на ряде модельных задач.
Задачи пороупругости для пластин рассмотрены во второй главе. Особенностью рассматриваемых математических моделей является то, что перемещения пористой среды описываются одной скалярной величиной — перемещением по нормали к срединной плоскости. В этих задачах пороупругости система состоит из эллиптического уравнения четвертого порядка для перемещений и параболического уравнения для давления. Отмечаются особенности конечно-элементной аппроксимации краевой задачи, строятся чисто неявные двухслойные схемы для системы уравнений пороупругости пластин. Основной результат состоит в построении схем расщепления по физическим процессам. Теоретическое рассмотрение дополняется результатами численных экспериментов.
В третьей главе разработанные схемы расщепления по физическим процессам применяются для численного решения задач пороупругости при разработке месторождений нефти и газа. В рамках приближения однофазной фильтрации рассматриваются задачи фильтрационной консолидации в трехмерном приближении. Вычислительный алгоритм базируется на конечно-элементной аппроксимации по пространству. Программное обеспечение для приближенного решения задач пороупругости разработано с использованием современных вычислительных технологий, свободных библиотек инженерных и научных вычислений, ориентированных на компьютеры параллельной архитектуры.
Численные эксперименты
Проблемам построения, обоснования и использования схем расщепления по физическим процессам (перемещение и давление) посвящено настоящее исследование. В первой главе обсуждаются проблемы численного решения краевых задач пороупругости (фильтрационной консолидации). Аналогичные математические модели используются при исследовании термоупругости. Базовая система уравнений включает стационарные уравнения Ламе для перемещений и нестационарные уравнения для давления или температуры. Вычислительный алгоритм основан на конечно-элементной аппроксимации по пространству. Для двухслойных схем с весами формулируются стандартные условия устойчивости. Вычислительная реализация таких схем основана на решении системы связанных уравнений для перемещений и давления (температуры). Строятся схемы расщепления по физическим процессам, когда переход на новый временной слой связывается с решением отдельных эллиптических задач для искомых перемещений и давления (температуры). Безусловно устойчивые аддитивные схемы строятся на основе выбора веса трехслойной схемы. Работоспособность предложенных схем расщепления продемонстрирована на ряде модельных задач.
Задачи пороупругости для пластин рассмотрены во второй главе. Особенностью рассматриваемых математических моделей является то, что перемещения пористой среды описываются одной скалярной величиной — перемещением по нормали к срединной плоскости. В этих задачах пороупругости система состоит из эллиптического уравнения четвертого порядка для перемещений и параболического уравнения для давления. Отмечаются особенности конечно-элементной аппроксимации краевой задачи, строятся чисто неявные двухслойные схемы для системы уравнений пороупругости пластин. Основной результат состоит в построении схем расщепления по физическим процессам. Теоретическое рассмотрение дополняется результатами численных экспериментов.
В третьей главе разработанные схемы расщепления по физическим процессам применяются для численного решения задач пороупругости при разработке месторождений нефти и газа. В рамках приближения однофазной фильтрации рассматриваются задачи фильтрационной консолидации в трехмерном приближении. Вычислительный алгоритм базируется на конечно-элементной аппроксимации по пространству. Программное обеспечение для приближенного решения задач пороупругости разработано с использованием современных вычислительных технологий, свободных библиотек инженерных и научных вычислений, ориентированных на компьютеры параллельной архитектуры.
Прикладное математическое моделирование пороупругости играет важную роль во многих областях науки и техники [16, 18, 105, 142]. В частности, необходимость учитывать деформационные процессы появляется в разработке нефтяных и газовых месторождений, гражданском строительстве, охране окружающей среды и биоинженерии.
Математические модели пороупругости включают уравнения упругости Ламе для перемещений среды и нестационарное параболическое уравнение для давления фильтрующейся жидкости [64, 142, 146]. Заметим, что математические модели термоупругости являются аналогичными[109, 118]. Наиболее важной особенностью математических моделей пороупругости является сильная связь между уравнениями системы. C одной стороны, уравнение для перемещений включает объемную силу, которая пропорциональна градиенту давления. С другой стороны, уравнение для давления включает сжимаемость среды, которая, в свою очередь, пропорциональна дивергенции скорости перемещений.
Вычислительные алгоритмы приближенного решения задач пороупруго-сти базируются, чаще всего, на конечно-элементной аппроксимации по пространству [35, 82, 97, 107, 108]. Для удовлетворения так называемого LBB-условия [50] мы не можем использовать конечно-элементные аппроксимации одного порядка для смещений и давления и должны ориентироваться на смешанные конечные элементы, на более низкий порядок аппроксимации для давления.
При аппроксимации по времени используются стандартные двухслойные схемы с весами. В работах [74] устойчивость и сходимость двухслойных схем для задач пороупругости при конечно-разностной аппроксимации по пространству исследуется в рамках общей теории устойчивости (корректности) операторно-разностных схем [22, 124].
Отдельного внимания заслуживают исследования по построению схем расщепления по физическим процессам, когда переход на новый временной слой осуществляется последовательным решением отдельных задач расчета перемещений и расчета давления. Различные классы таких схем строятся на основе аддитивного расщепления оператора задачи [103, 138]. Для задач пороупругости схемы расщепления строятся и используются в работах [36, 89–91, 106].
Несмотря на длительную историю и большое внимание, в настоящее время мы не имеем строгих математических результатов об устойчивости схем расщепления по физическим процессам для задач пороупругости и термоупругости. До сих пор выводы об устойчивости и работоспособности таких схем делаются не на основе их теоретического исследования, а на базе анализа вычислительных экспериментов по приближенному решению некоторых модельных задач. В нашей работе дается попытка восполнить этот пробел.
В работе на основе принципа регуляризации операторно-разностных схем А.А. Самарского [21, 22] строятся схемы расщепления по физическим процессам для задач пороупругости и термоупругости. Успех достигается за счет перехода к трехслойной разностной схеме, когда часть членов системы уравнений берется с нижнего временного слоя и выбора весового параметра (параметра регуляризации). 1.2. Постановка задачи
Рассмотрим задачу пороупругости в ограниченной области c границей . Математическая модель содержит уравнения для перемещения и давления, которые получаются из уравнения равновесия и неразрывности, соответственно. Обозначим вектор смещения пористой среды через и , а давление фильтрующейся жидкости через р. В квази-стационарном приближении напряженно-деформированное состояние пористой среды описывается с помощью уравнения равновесия [42]:
Численное решение задачи прогиба пластин
Рассматриваются проблемы численного решения краевых задач поро-упругости (фильтрационной консолидации) для пластин. Базовая система уравнений включает двумерное бигармоническое уравнение для вертикальных перемещений и нестационарное уравнение для давления в пористой среде. Вычислительный алгоритм основан на конечно-элементной аппроксимации по продольным координатам и разностной аппроксимации по времени. Для двухслойных схем с весами формулируются стандартные условия устойчивости. Вычислительная реализация таких схем основана на решении системы связанных уравнений: эллиптического уравнения четвертого порядка для перемещений и эллиптического уравнения второго порядка для давления. Строятся безусловно устойчивые схемы расщепления по физическим процессам, когда переход на новый временной слой связывается с решением отдельных задач для искомых перемещений и давления.
При прикладном математическом моделировании большое внимание уделяется задачам пороупругости, задачам фильтрационной консолидации [105, 142]. Базовые математические модели включают уравнения Ламе для перемещений среды и нестационарное параболическое уравнение для давления фильтрующейся жидкости. Подобные математические модели используются при описании термоупругих напряжений [109, 118], когда стационарные урав нения для перемещений дополняются нестационарным уравнением теплопроводности. Наиболее важная особенность математических моделей пороупру-гости и термоупругости является того, что уравнения системы завязаны между собой. В задачах пороупругости (термоупругости), с одной стороны, уравнение для перемещений включает объемную силу, которая пропорциональна градиенту давления (температуры). С другой стороны, уравнение для давления (температуры) включает член, который описывает сжимаемость среды (дивергенция скорости перемещений).
Во многих прикладных задачах имеет место малость поперечного размера расчетной области по отношению к продольному размеру. В этом случае используются упрощенные модели: математические модели упругих пластин и оболочек [57, 120]. Понижение размерности задачи приводит к тому, что перемещение пластины в поперечном направлении описывается эллиптическим уравнением четвертого порядка. Учет неоднородности давления (температуры) по пластине в задачах пороупругости (термоупругости) проводиться в различных приближениях. Часто можно не учитывать неоднородность по толщине и описывать давление (температуру) в рамках двумерного приближения [52, 96].
Более содержательные модели пороупругости и термоупругости пластин и оболочек базируются на учете неоднородности давления и температуры по толщине пластины [130, 134]. В этом случае давление (температура) по всей пластине (трехмерная задача) определяется из уравнения фильтрации. Для нестационарных задач термоупругости характерно зацепление уравнений: ис-точниковый член в уравнениях упругости зависит от давления, а в уравнении фильтрации имеется член, зависящий от динамики перемещений упругой среды.
При аппроксимации по времени для задач пороупругости и термоупру-59 гости используются стандартные двухслойные схемы с весами. В работах [74, 99] устойчивость и сходимость двухслойных схем для задач пороупру-гости при конечно-разностной аппроксимации по пространству исследуется в рамках общей теории устойчивости (корректности) операторно-разностных схем [22, 124].
Отдельного внимания заслуживают исследования по построению схем расщепления по физическим процессам, когда переход на новый временной слой осуществляется последовательным решением отдельных задач расчета перемещений и расчета давления (температуры). Различные классы таких схем строятся на основе аддитивного расщепления оператора задачи [103, 138]. Для задач пороупругости схемы расщепления строятся и используются в работах [36, 89, 90, 94, 106].
В своем исследовании мы рассматриваем проблемы численного решения задач пороупругости (термоупругости) для пластин. Используется конечно-элементная аппроксимация по продольным координатам. Отдельно обсуждаются проблемы численного решения краевых задач для бигармонического уравнения. Исследуется устойчивость двухслойных схем аппроксимации по времени. Основное внимание уделяется построению схем расщепления по физическим процессам для задач пороупругости и термоупругости для пластин. Успех достигается за счет перехода к трехслойной разностной схеме, когда часть членов системы уравнений берется с нижнего временного слоя.
— это тело, имеющее форму прямой призмы с небольшой по сравнению с размерами оснований толщиной h. Обозначим через Q — поперечное сечение пластины и пусть жі, #2 — коорди наты вдоль пластины, х% — координата поперек пластины, —/г/2 х% /г/2.
Уравнения пороупругости (1.9), (1.10) описывает общий случай, когда мы имеем дело с трехмерным телом. Но процесс деформирования пластин обладает некоторыми особенностями [57]: нормаль к срединной плоскости (плоскость, делящая пластину пополам) при изгибе не искривляется и не меняет своей величины, а поворачивается, оставаясь ортогональной к срединной плоскости, нормальные напряжения между слоями пластины, параллельными срединной плоскости, пренебрежимо малы по сравнению с нормальными напряжениями в сечениях, ортогональных к этой плоскости, при изгибе пластины срединная плоскость переходит в искривленную поверхность, не испытывая осевых и сдвиговых деформаций, прогиб срединной плоскости мал относительно толщины пластины.
Далее на основе данных допущений получим систему уравнений, описывающую процесс деформирования пористой насыщенной жидкостью пластины. Из вышеупомянутых допущений следует, что, во-первых, компонентами тензоров деформаций и напряжений єіз, 23 и с"зз пренебрегаем и, во-вторых,
Вычислительный алгоритм
Расчеты проводились на трех вычислительных тетраэдальных сетках с локальным сгущением вокруг скважины, которая вырезана из расчетной области (рис. 31). Сетки построены с помощью программы Gmsh. Количество вершин, ячеек и количество неизвестных в полях давления и смещения для каждой сетки представлены в табл. 13.
Для решения линейной системы алгебраических уравнений используем обобщенный метод минимальный невязок (GMRES) с неполной ILU факторизацией в качестве предобуславливателя. Расчеты проводились на вычислительном кластере Северо-Восточного федерального университета имени М.К. Аммосова Ариан Кузьмин.
В табл. 14 представлена зависимость общего времени расчетов в секундах от количества запущенных процессоров и сетки при использовании неявной схемы (3.8), (3.9) и схемы расщепления (3.6), (3.7). Здесь, в общее время входят затраты на инициализацию расчетов, загрузку и разделения расчетной сетки, подготовку и сборку матрицы и вектора линейной системы, а также расчет 10 временных шагов при шаге по времени = 1 день. В табл. 15 -17 представлены время счета и количество итераций для решения линейной системы связанной задачи, задачи для давления и смещения, соответственно.
Видно, что схема расщепления значительно сокращает время расчета (около 60%). Также, при согласовании размера используемой сетки и числа процессоров обеспечивается достаточно высокая степень эффективности параллелизации. При относительно малом количестве ячеек использование большого количества процессоров может замедлить процесс расчета, так как время расходуется на обмен информацией между процессорами.
Далее оценим величину параметра и веса , при которой схема расщепления (3.6), (3.7) безусловно устойчива согласно теореме 3, глава 1. Для этого решим спектральную задачу (1.76). Для рассматриваемых входных параметров (табл. 12) получим следующие результаты:
Для исследования устойчивости схемы расщепления будем смотреть значение давления вблизи скважины в точке (450.0, 500.0, 200.0), чтобы уловить начальные колебания давления в случае неустойчивости. На рис. 32 представлены графики динамики давления при различных значениях веса для шагов по времени = 1 и 2 дня. Стандартная схема расщепления с = 1.0 (схема с дренированием) неустойчива для любого шага по времени. Заметим, что величина шага по времени не влияет на устойчивость. При увеличении веса до 1.4 схема становится устойчивой, что соответствует теоретическим результатам. Однако, в начальный момент времени все еще наблюдаются колебания значения давления. Они полностью исчезают при = 1.8.
Далее, для исследования зависимости точности решения от вычислительных параметров (сетки и шага по времени) будем смотреть величину пластового давления на краю рассматриваемой области в точке (0.0, 500.0, 200.0). На рис. 33 представлен график динамики давления для различных сеток. Расчеты проводились при = 1 день и = 1.8. В начальный момент времени качество сетки не влияет на значение давления. К концу расчета (100 дней) давления для разных сеток немного отличаются. Чем подробнее сетка, тем больше значение давления. Аналогичные результаты получены для зависимости давления от шага по времени (рис. 34). Видим, что чем меньше шаг по времени, тем больше давление.
Наконец, на рис. 35 и 36 показаны распределения давления и вертикальной компоненты вектора смещения в различные моменты времени (1 месяц, 1 год, 2 года, 4 года, 6 лет), соответственно. Данное решение получено с помощью схемы расщепления с весом = 1.8 на сетке 3 при = 1 день. Смещение представлено на деформированной области и умножено 100 для наглядности. Расчеты проводились на 64 процессорах. Время расчета составило 88 мин. Декомпозиция расчетной сетки представлена на рис. 37. За 6 лет разработки пластовое давление уменьшилось до 34.3 МПа и кровля пласта опустилась на почти 1 м в области скважины и на 0.93 м на границах.
Далее рассмотрим случай, когда для поддержание пластового давления и предотвращения опускания поверхности пласта применяется площадное заводнение. Для это в дополнению к добывающей скважине пробурены четыре нагнетательные скважины (рис. 38). Данные скважины имеют радиус г = 0.1 м и вскрывают пласт не на всю толщину а на h = 150.0 м от подошвы пласта. Закачка производится при забойном давлении, которое линейно возрастает от Pi = 40.0 до 45.0 МПа. Все остальные параметры остаются такими же как в предыдущем случае (табл. 12).
Проведем численное моделирование процесса пороупругости, используя схему расщепления (3.6), (3.7) с весом в = 1.8 и временным шагом г = 1 день на вычислительной сетке, содержащей 484948 ячеек (рис. 39). Результаты расчетов показаны на рис. 40 и 41. Смещение представлено на деформированной области и умножено 100. За два года добычи поверхность пласта около скважины опустилась на 0.4 м, но в результате заводненения пластовое давление увеличилось и в конце расчетного периода (6 лет) уровень опускания поверхности уменьшился до 0.09 м. В то же время, в районе нагнетательных скважин поверхность приподнялась на 0.32 м.
