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

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

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

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

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

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

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

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

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

1.5.1. ОСНОВНЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ

1.5.1.1. Рассмотрим следующее дифференциальное уравнение параболического типа

(1.5.1)

Здесь под x, y- подразумеваются пространственные координаты, t-время, U -искомая функция (P или P ² в случае газа); s, g -коэффициенты уравнения, зависящие от координат x, y времени t и самой функции U (в линейном случае постоянные), q-мощность источника (стока), приходящейся на объем ячейки.

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

(1.5.2)

(i, j=1,2…M, N)

где

Отметим, что рассматривается пока прямоугольная сеточная область. Легко видеть, что система (1.5.2) образует пяти-диагональную матрицу.

H D B E F i +1 i- 1 j i j +1 j -1

Рис.1.2

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

(метод Гаусса, Жордана, Холецкого и пр.) для системы (1.5.2) непригодны, поскольку требуют порядка 0 операций (h -шаг сетки). Использование формул Крамера также исключается, т.к. для вычисления определителя N -го порядка требуется ~ N! умножений.

Значительно более простыми в отношении программирования являются итерационные

методы. и, несмотря на то, что они не приводят за конечное число шагов к точному реше-

нию, являются наиболее предпочтительными.

Однако простейшие итерационные методы (Гаусса-Зейделя, Якоби) сходятся сравнительно медленно (~ операций или 0()).

В 1950 Франкел и Янг (независимо) предложили специальный итерационный метод – “экстраполяции Либмана” или “метод последовательной верхней релаксации” (Янг), которые требовали ~ 0() итераций. Затем появилось много работ по итерационным методам решения специальных линейных уравнений типа (1.5.2) с применением полиномов Чебышева, которые требуют также ~ 0() итераций. Наконец в 1955г. Шелдон предложил метод симметричной последовательной верхней релаксации, в котором сочетаются метод последовательной верхней релаксации и метод с использованием полиномов Чебышева, который требовал итераций. В том же 1955г. появился метод переменных направлений (Дуглас и Писмен и Речфорд), требующий итераций.

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

Есть еще одна группа итерационных методов – вариационные методы: метод наискорейшего спуска, метод сопряженных градиентов. Эти методы построены на принципе минимизации соответствующей квадратичной формы. (См. например работу Саульева [25], или

Самарского А.А. [26]).

Мы не будем рассматривать все эти итерационные методы, а остановимся более подробно на методе переменных направлений и на его дальнейшей модификации - суммарной аппроксимации (ЛОС) Самарского А.А.

1.5.1.2. Прежде чем перейти к методу переменных направлений рассмотрим метод матричной прогонки [14, 25].

Рассмотрим вновь уравнение (1.5.2), которое применено к прямоугольной области (M*N) На границе задается равенство функции нулю, т.е.

(1.5.3)

Тогда в матричной форме можно записать(1.5.2) так,

(1.5.4), где

; ;

Пусть решение ищется в виде:

(1.5.5)

По аналогии с обычной прогонкой имеем:

, (1.5.6)

, (1.5.7)

(i= M-1, …, 1)

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

1.5.2.МЕТОД ПЕРЕМЕННЫХ НАПРАВЛЕНИЙ И ЛОКАЛЬНО-ОДНОМЕРНАЯ

СХЕМА.

Как уже указывалось выше, экономичными разностными методами называются такие, у которых число итераций не превышает 0(). Все экономичные методы имеют одну общую алгоритмическую идею: процесс отыскания приближенного решения многомерной задачи разбивается на несколько этапов, на каждом из которых решается простая задача. Так, для уравнения второго порядка такой простой, ² первичной ² задачей является трех точечная разностная задача, которая решается методом прогонки. Эта трех точечная задача, как правило, может быть трактована как разностная аппроксимация одномерного (по X_a) дифференциального уравнения. Таким образом, экономичный алгоритм решения сложных задач есть цепочка простых алгоритмов.

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

1.5.2.1. Основные черты метода.

1.Переход со слоя (n) на слой (n+1) осуществляется при помощи последовательности обычных (двухслойных, трехслойных) схем.

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

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

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

Вместо разностного аналога уравнения (1.5.1) по формуле (1.5.2) будем иметь следующую систему разностных уравнений:

(1.5.8a)     (1.5.8б)

Здесь наряду с основными значениями искомой функции вводится промежуточное значение, , которое формально рассматривается на момент времени . Переход от временного слоя n к слою n+1 совершается в два этапа с шагами 0.5t.

При i= 0,M и j= 0,N задаются граничные условия, например непроницаемая граница.

Отборы из скважин могут быть заданы в виде источников .

Система (1.5.8) имеет (N+1)(M+1) неизвестных и состоит из (M-1)(N-1) уравнений, к которым для замыкания системы следует добавить (2M+2N) условий на сеточной границе.

Характерным для (1.5.8а) и (1.5.8б) является, то, что они имеют трех диагональную матрицу и могут быть решены методами прогонки.

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

Итак, методом переменных направлений сначала решается система (1.5.8a). Прогонка осуществляется вдоль оси X. При этом значения второй производной по Y берется с нижнего временного слоя (n). Решение осуществляется для каждой строки. Затем решается система (1.5.8б) для каждого столбца вдоль оси Y.

Получаемое таким образом решение является искомым и соответствует (n+1) временному слою. Таким же образом отыскивается решение задачи на каждом следующем временном слое.

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

1.5.2.2. Сейчас, мы рассмотрим разностную схему, которую Самарский А.А. называет локально-одномерной схемой (ЛОС) [26]. Эта схема легко обобщается на случай переменных и разрывных коэффициетов и для произвольной области, а также на многомерный случай.

Для этой схемы разностные уравнения для дифференциального уравнения (1.5.1) принимают вид:

(1.5.9)

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

Однако здесь

(1.5.10)

т.е. сумма значений мощности источника на n+1/2 и n+1 временных слоях равна искомому значению мощности источника.

Необходимо отметить, что при использовании ЛОС имеет значение процесс организации прогонки и момент вывода результатов на печать. Решение систем разностных уравнений начинается одной серией прогонок, на пример по оси X. За тем выполняются дважды прогонки по оси Y, дважды по оси X и т. д. Перед выводом на печать результата производится лишь одна серия прогонок по оси X. При выполнении этого правила расчет дает правильные не только количественные, но и качественные результаты. Если это правило не выполняется, то при верных, например, забойных давлениях получаются некачественные карты изобар, построенные по совокупности всех узловых точек области интегрирования.

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

1.5.3 ЗАДАНИЕ ГРАНИЧНЫХ УСЛОВИЙ ДЛЯ ДВУХМЕРНОЙ ЗАДАЧИ.

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

1.Внешний контур проводится между соседними узлами сетки.

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

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

Рис.1.3

; ;

(1.5.11a)

Т.е. () и () в уравнении (1.5.2)

На скважинах задаются условия, что мощность источников ()- функция времени, т.е.

(1.5.11б)

Т.к. мощность источника, это дебит, приходящийся на объем элемента пласта с размерами

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

Действительно, пусть U _ij = , тогда (1.5.12)

Где или , при этом ,