Уравнения типа уравнения теплопроводности и фильтрации

Рассмотрим представление в конечно – разностной форме одномерного （一维）дифференциального уравнения параболического（抛物线）типа

; ,

(1.2.1)

которое описывает фильтрацию сжимаемой（压缩）жидкости (аналогично уравнению теплопроводности热传导性)

Используя уравнения (1.1.5) и (1.1.7) предыдущего раздела имеем выражение

(1.2.2)

Здесь - погрешность（误差，修正量） аппроксимации исходного дифференциального уравнения (1.2.1) конечно – разностным уравнением. Принимается, что

В уравнении (1.2.2) левую часть можно рассматривать на различных временных уровнях: либо на слое j Dt и тогда имеем уравнение

(1.2.3)

либо на временном слое (j+1) Dt и тогда уравнение (1.2.2) имеет вид

(1.2.4)

При записи этих уравнений пренебрегается величиной

В уравнении (1.2.3) имеется лишь одна неизвестная величина (Считается, что все величины на временном слое известны). Такое уравнение называется явным сеточным уравнением.

Уравнение (1.2.4), где имеются три неизвестные величины , , называется неявным.

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

Записывая неявное уравнение (1.2.4) для точек , получаем систему из (n-1) уравнения с (n+1) неизвестным. Граничные условия в точках i=0 и i=n дают еще два условия. Следовательно, чтобы решить задачу на временном слое , требуется решить систему (n+1) уравнения с (n+1) неизвестным.

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

1.2.2. Использование сеточных методов для решения краевых（边缘的，低于最大限度的） задач приводит к рассмотрению вопросов сходимости收敛性 и устойчивости 稳定性 их.

Необходимость рассмотрения сходимости и устойчивости сеточных методов связана:

1) с погрешностью误差 аппроксимации近似值 дифференциальных уравнений и соответствующих相关的краевых 边缘условий конечно-разностными уравнениями;

2) с погрешностью решения на каждом временном слое .

Если сеточный метод, дает такое решение, которое при изменении шагов ∆x и ∆t (при , ) (причем , т.е. пространственные（空间的） приращения（增量）стремятся к нулю по определенной зависимости от временного приращения ∆t, когда последнее стремится 趋向к нулю), стремится к точному решению задачи, то такой метод называется сходящимся, а конечно-разностное уравнение согласуется с соответствующими дифференциальным уравнением в частных производных.

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

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

Если рассматривать разностные методы решения дифференциальных уравнений с точки зрения линейных операторов в векторном пространстве, то всякое решение на временном слое, можно выразить через решение на временном слое, применив к последнему, оператор（算子）, называемый матрицей (множителем) перехода от одного временного слоя к другому. При этом устойчивость（稳定性） требует, чтобы, матрица（矩阵） перехода была равномерно ограниченной.

Необходимым условием устойчивости по Нейману является условие, что , при ; и всех , где – сетка пространства, – радиус спектра матрицы [22].

Пусть , тогда , и в частности . Для из интервала величина ограничена линейным выражением . Отсюда получаем условие:

max , (a)

при , -порядок (число всех слоев по и всех ).

При этом , а - собственное значение матрицы перехода .

В простейшей задаче теплопроводности (1.2.1) условия устойчивости имеет вид

(б)

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

1.2.3.2. Собственное значение матрицы перехода связано с соотношением величин шагов и

В случае явного уравнения (1.2.3) решение является устойчивым, если существует следующее соотношение между пространственными и временными шагами:

(1.2.5)

В случае неявного разностного уравнения (1.2.4) никаких ограничений на величины шагов и не накладываются, т. е. говорят, что уравнение (1.2.4) абсолютно устойчиво.

1.2.4. Однако это не означает, что при пользовании неявным методом допустимы любые шаги на осях Dx и Dt.

Здесь мы сталкиваемся с понятием ошибки аппроксимации.

Действительно, вычтем из (1.2.1) (1.2.2) и получим после преобразования выражение:

(1.2.6)

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

1.2.5. Если является точным решением дифференциального уравнения, то величина, стоящая в скобках в уравнении (1.2.6), обращается в нуль. Кроме того, членами порядка Dt и (Dx)² являются, согласно формуле Тейлора

(1.2.7)

Т.к. величина давления удовлетворяет уравнению (1.2.1), то она также удовлетворяет и дифференциальному уравнению:

(1.2.8)

и главным членом в правой части (1.2.6) будет:

Если, при этом, Dt и Dx связаны между собой соотношением

, (1.2.9)

то главный член правой части уравнения (1.2.6) обращается в нуль. Погрешность аппроксимации в этом случае становится и

1.2.6. Использование явного сеточного уравнения (1.2.3) возможно в том случае, если соблюдается условие (1.2.5). Это ограничение на шаг по времени оказывается очень жёстким. Для соблюдения условия устойчивости шаг Dt приходится брать очень малым. Это увеличивает общее число шагов по времени, а, следовательно, и общий объём вычислительных работ. В связи с этим, несмотря на достаточную простоту явного сеточного метода, его использование на практике весьма ограничено.

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

Выбор шага по пространственной переменной на практике часто осуществляется экспериментальным путём. Для этого на ЭВМ проводятся вычисления с некоторыми шагами D x. Затем вычисления проводятся с шагом D x /2. Если оказывается, что решения отличаются на величину меньшую заданной погрешности e, то шаг D x считается незавышенным для достижения заданной точности. В противном случае делается просчёт с меньшим шагом, D x /4 и т.д.

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

Например, замена обычных координат в осе симметричном случае на логарифмические.

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

Один из алгоритмов увеличения шага по времени таков. Рассчитываются с малым начальным шагом D t по времени 2 временных шага. Затем с шагом просчитывается один шаг по времени. Полученные два решения на момент времени сопоставляются. Если эти решения различаются на величину большую, чем заданная погрешность , то дальнейший счет ведется с шагом . В противном случае расчет ведется с шагом и т.д.

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

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

1.2.7. В заключении приведем ряд конечно-разностных схем применительно к уравнению (1.2.1)

, s = const > 0

Таблица 1

	Мнемонические схемы	Разностные уравнения и погрешность аппроксимации
	1		Явная схема устойчива, если
	2		Неявная схема всегда устойчива
	3		Частный случай 1
	4		Неявная схема всегда устойчива
	5		при устойчива, если , при всегда устойчива.
	6		Разностное уравнение такое же как в случае 5,но   , схема устойчива.
	7		- схема всегда неустойчива.
	8		Схема всегда устойчива.
	9		Схема всегда устойчива.
	10		Схема всегда устойчива, включает случаи 3 и 9.
	11		Разностное уравнение как в случае 10, но Схема всегда устойчива
	12		Схема всегда устойчива
	13
	14		Если , то имеем обобщение схемы 12 (например) Схема устойчива.

1.3. ПРОГОНКА.

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

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

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

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

Это по существу специальный метод исключения Гаусса. Идея заключается в замене разностного уравнения второго порядка тремя разностными уравнениями первого порядка.

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

(1.3.1)

Начальные и граничные условия:

t = 0 P= P_k = const (1.3.2)

r= R₀ (1.3.3a)

r= R _k (1.3.3б)

Далее пусть уравнение состояния для газа имеет вид:

, (1.3.4a)

Если газ идеальный, то имеем ,

если газ реальный, то имеем ,

для упругого пласта имеем (1.3.4б)

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

;

(1.3.5)

С учетом (1.3.5) уравнение (1.3.1) и граничные условия (1.3.2) и (1.3.3) имеют вид (для простоты в дальнейшем звездочки будем опускать.)

	(1.3.1’)
	(1.3.2’)
	(1.3.3’a)
	(1.3.3’б)

Уравнение (1.3.1) является общим нелинейным параболическим дифференциальным

уравнением. Из него легко получить, например, дифференциальные уравнения для течения идеального газа. Для этого надо положить , а =1.

1.3.3. Представим теперь задачу (1.3.1’)-(1.3.3’б) в конечно-разностной форме.

При этом рассмотрим пространственно-временную непрерывную область в виде сеточной области.

;

¯ ¤

J+1 J j-1									Будем использовать неявную разностную схему вида крест.   Граница пласта и скважина расположены в середине интервалов, чтобы получить
0 1 2 I- 1 I I+ 1 m- 1 m Рис.1.1

лучшую аппроксимацию, порядка

	(1.3.6)
	(1.3.7)

Для записи разностной схемы уравнения (1.3.1’) поступим следующим образом. Обозначим тогда левая часть уравнения (1.3.1^’) записывается в виде и далее, заменяя на разности , получим

Правую часть (1.3.1^’) представим согласно мнемонической схемы на рис.1.1 и, учитывая, что пространственные разности берутся на верхнем временном слое, имеем после некоторых преобразований

Или окончательно

	(1.3.8)
(i =1, …n-1)

Уравнение (1.3.8) совместно с (1.3.6) и (1.3.7) представляют систему (n+1) уравнения с (n+1) неизвестным. При этом если коэффициенты при неизвестных будут постоянными (а это на каждом временном слое или итерации можно получить, определяя коэффициенты и с предыдущего слоя или итерации), то система линейна и имеет три диагонали.

1.3.3.1.Перепишем в общем случае задачу (1.3.8), (1.3.6), (1.3.7) так

	(1.3.9)
;	(1.3.10)

- общий вид граничных условий.

Пусть имеет место рекуррентное соотношение:

,

(1.3.11)

Подставим в (1.3.9) выражение , тогда

(1.3.12)

Сопоставляя (1.3.11) и (1.3.12), видим, что они совпадают при любых значениях , если:

;

(1.3.13)

Если коэффициенты в точке известны, то при помощи формул (1.3.13) можно просчитать коэффициенты , и т. д.

Затем, зная значения по формуле (1.3.11), можно определить значения для всех .

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

Тогда будем иметь

(1.3.14)

В нашем конкретном случае имеем (см.1.3.6.)

Итак, для функций E_i и F _i имеем формулы прямой прогонки (1.3.13), (1.3.14)

Если E _n-1 и F _n-1 известны, то из граничного условия P _n = ₂ P _n-1+ µ ₂

и P _{n- 1} = E_n-1 P _n-1+ F _n-1 можно получить значение P_n на границе. Т.е.

(1.3.15)

В нашем случае с учётом (1.3.7) имеем

Итак, формулы (1.3.11) и (1.3.15) позволяют определить значения P _i для всех значений

i = n, n-1…0. (Это формулы обратной прогонки.)

1.3.3.2. Устойчивость метода прогонки.

Формулы прогонки (1.3.13), (1.3.14) и (1.3.15), (1.3.11) были нами получены формально. Мы делили на выражения (b _i- c _i E_{i- 1}) и (1- ₂E_n-1), предполагая, что они не равны нулю. Сейчас мы укажем достаточные условия, когда это можно делать.

Пусть

(1.3.16)

Для устойчивости метода прогонки достаточно иметь для всех i=1,2,..n-1

Это действительно так.

Рассмотрим разность

, при .

Поскольку , то , т.е.

.

Отсюда, видно, что , если . Тогда все при .

Рассмотрим теперь неравенство

, т.к. или , или ,

т.е.

Таким образом, при выполнении условий (1.3.16) разностная задача, решаемая методом прогонки, имеет единственное решение.

1.3.4. В силу того, что решение ведется на ЭВМ приближенно, с конечным числом значащих цифр, возникают ошибки округления. Из-за них фактически находится не функция - решение задачи (1.3.9), (1.3.10), а - решение той же задачи с возмущенными коэффициентами и правыми частями . При этом, если процесс вычислений происходит с возрастанием ошибки округления, то это может привести как к потере точности, так и к невозможности продолжить вычисления из-за роста получаемых величин.

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

, где - ошибка округления.

Для погрешности имеем уравнение

,

или

,

откуда видно, что ошибка при экспоненциально возрастает с ростом i.

Для метода прогонки, в силу условия , ошибка при определении не нарастает.

Действительно, из уравнений

следует

, , т.е.

, т.к.

Отсюда и следует устойчивость метода прогонки.

Замечание. Если учесть, что в ходе вычислений возмущающимися являются и коэффициенты и , то показывается, что ошибка в определении в задаче (1.3.9),(1.3.10) пропорциональна квадрату числа узлов

,

где - ошибка округления.

Отсюда видно, что существует связь между точностью определения решения задачи, числом уравнений и числом значащих цифр ЭВМ, так как

(Если разрядность машины 12 значащих цифр - , n = 100, тогда точность решения )

1.3.5. Выше уже отмечалось, что уравнение (1.3.8) описывает нелинейную фильтрацию газа в круговом пласте. Для решения задачи численным методом необходимо линеаризовать нелинейные коэффициенты. Этого можно достичь с использованием их значений, вычисленных при значении давлений на предыдущем временном слое, или для более точного учета их необходимо итерировать на каждом временном слое. Первым приближением служат нелинейные коэффициенты, вычисленные по данным решения на j -ом временном слое. Затем определяется приближенное решение на (j +1)-ом временном слое. Во втором приближении для вычисления нелинейных членов используется решение на (j+1)-ом слое. Вновь определяется решение на (j +1)-ом временном слое и т. д., до тех пор, пока не будет выполнено неравенство

Очевидно, что, начиная счет в обратном порядке, получаем левуюпрогонку. В этом легко убедиться самостоятельно. Если у нас имеется какой-либо внутренний узел , то можно получить метод встречных прогонок, комбинируя правую и левую прогонки. Отметим, что все эти прогонки являются устойчивыми. (См. Самарский А. А. ‘’Теория разностных схем’’. Наука, 1972 [26]).

1.4. ПОТОКОВАЯ ПРОГОНКА.

Этот вариант прогонки был разработан для решения задач с сильноменяющимися от точки к точке коэффициентами.

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

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

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

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

Рассмотрим краевую задачу (1.3.5), (1.3.6), (1.3.7), которая равносильна разностной краевой задаче (1.3.9), (1.3.10) предыдущей части, т.е.

(1.4.1) (1.4.2)

Где

(1.4.3)

Формулы обычной прогонки для (1.4.1), (1.4.2) с учетом (1.4.3) принимают вид:

, (i =0,1,2…n-1)

(1.4.4)

Введем новую неизвестную разностную функцию (поток) по формуле

(1.4.5)

И перепишем (1.4.1) и краевые условия (1.4.2) в виде

(1.4.6)   (1.4.7а) (1.4.7б)

Подставляя значение из (1.4.5) в первое уравнение (1.4.4) получим

(1.4.8)

Вводя обозначения ;

перепишем (1.4.8) следующим образом

(1.4.9)

Исключая из (1.4.6) и (1.4.9), получим

(1.4.10)

(i =n-1,…,2,1)

Для определения и напишем рекуррентные соотношения

, или

(1.4.11)

(i=1,2,…,n-1)

или

(1.4.12)

(i =1,2,…,n-1)

При соблюдении условий (1.4.3) из (1.4.11) следует, что . Тогда коэффициент в формуле (1.4.10) всегда меньше единицы, что обеспечивает устойчивость при вычислении .

Из сравнения первого краевого условия (1.4.7a) c (1.4.9) при I=1 находим

(1.4.13)

Для определения P воспользуемся следующими формулами

a

(1.4.14)

(1.4.15)

Из (1.4.14) и (1.4.15) видно, что прогонка устойчива. Для счета по формулам (1.4.10), (1.4.14) и (1.4.15) нужно знать величины и P _n. Они могут быть определены из второго краевого условия (1.4.7б) и соотношения (1.4.9) при i =n, при a <1

P

(1.4.16)

и при a

P

(1.4.17)

Из (1.4.3) следует, что знаменатели формул (1.4,16), (1.4.16а) всегда больше нуля.

Итак, алгоритм потоковой прогонки следующий:

1.)Вычисляем и по формулам (1.4.13)

2.) Для i =1,2… n -1 последовательно находим β_{i +1} по формуле (1.4.11) (первая часть a_i < 1, вторая часть a_i ³ 1) и g_{i+ 1} по формуле (1.4.12) (первая часть a_i < 1, вторая часть a _i ³ 1)

3.) Определяем P_n и W _n по формулам (1.4.16), если a _i < 1 и по формулам (1.4.16a), если a_{n- 1} ³1

4) Для i = n -1,…1, вычисляем W_i по формуле (1.4.10) и P _i по формуле (1.4.14) при a _i ³1 и по формуле (1.4.15) при a_i < 1

Замечание. Выше были приведены формулы для вычисления не только функции P _i, но и потока W _i. При больших значениях a _i вычисление потока по формуле W_i = a_{i -1}(P_{i -1} - P_i) приводит к существенной потере точности. Это и послужило одной из причин введения потока W _i в качестве дополнительно искомой функции и вычисления её по формуле (1.4.10).