Численное моделирование и анализ результатов расчета задачи фильтрации газа в пористой среде



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

Ключевые слова: фильтрация, пористая среда, метод прогонки, квазилейная, давления, пласт, вязкость.

Рассматривается краевая задача фильтрации газа в пористой среде, математическая модель которой описывается нелинейным дифференциальным уравнением к параболического типа второго порядка с переменнымы коэффициентами [1]:

(1)

Уравнения (1) интегрируется при следующих начальных и граничных условиях:

p(x,y,t0)=pH (х,y),(x,y)G;(2)

;(3)

.(4)

Здесьn – внутренняя нормаль к границе Г, qi – дебит i-й газовой скважины,

si– контур i-й скважин, p – давление в области газоносности, pH – начальное давление в области газоносности, k– коэффициент проводимости пласта, m– коэффициент пористости пласта, – коэффициент вязкости газа, n1– нормаль к контуру s.

Обезразмерывание переменнных осуществляются по следующим формулам:

p= p0p*; k=k0 k*; μ= μ0 μ*; x=L x*; y=L y*;

.

где приняты следующие обозначения: p0 – характерное значение давления в области газоносности, k0– характерное значение проницаемости пласта, μ0– характерное значение вязкости газа, L – характерная длина пласта.

Тогда задачи (1) – (4), можно переписать в следующем виде:

(5)

p(x,y,t)=pH (х,y),t=t0(x,y)G(6)

(7)

(8)

Безразмерное дифференциальное уравнение (5) с краевыми условиями (6)-(8) решаются численно. При этом для получения разностной задачи используется алгоритмическая идея неявной схемы переменных направлений [2] и алгоритм метода прогонки по каждому направлению.

Для перехода от дифференциальной задачи фильтрации к разностной, области G+Г разбивается на элементарные ячейки  строил которых параллелные осям OX и OY.

гдеNj – число узлов на прямой yj, Mi–число узлов на прямой xi, – шаг сетки.

Тогда, если использовать схему переменных направлений по переменной х для дифференциального уравнения в частных производных (5), описывающую процесс фильтрации, в области можно заменить следующим разностным оператором при +0,5:

(9)

где – значения функции давления k временном слое;

– значения функции давления k+0.5 временном слое;

Очевидно, что полученные разностные уравнения относительно функции давления pнелинейные. Здесь для решения разностной задачи применяется итерационный метод, основанный на методике квазилинеаризации нелинейных членов [3]. Согласно этому методу нелинейные члены разностного уравнения (9) представляется в виде:

. (10)

Здесь,– приближенное значение функции p, которое уточняется в процессе итерации , при этом .

Итерационной процесс продолжается до тех пор, пока не выполнится условие

(11)

где  – точность итерации, заранее известная, заданная малая величина, s – номер итерации.

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

.

Тогда вместо разностных уравнений (9) можно записать следующие квазилинейные разностные формулы для k+0.5 временного слоя

Отсюда можно записать следующие трёхточечные разноcтные уравнения

.(12)

где

Теперь используя схему переменных направлений по переменной y для дифференциального уравнения в частных производных (5), описывающего процесс фильтрации в области, можно заменить следующим квазиленейными разностными операторами, при k+1 – временном слое:

где – значения функции давления k+1 временном слое;

– значения функции давления k+0.5 временном слое;

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

(13)

где

Для численной модели краевой задачи (1)–(4) разработана программное обеспечения (ПО) по алгоритму из следующих этапов.

На первом этапе вычислений основных показателей разработки газового месторождения осуществляется вычислением значений функции давления k+0.5 временном слое по направлению переменной x при фиксированной переменной y. При этом применяется метод прогонки по направлению переменной x, вычисление осуществляется следующим образом:

-вычисления коэффициентов ai, bi, ci, di(i=1,2..Nj) трёхточечного разноcтного уравнения (12).

-определение первого значения прогоночных коэффициентов 0и 0 из граничных условий левой части дискретной области фильтрации.

-вычисления значения прогоночных коэффициентов iиi (i=1,2..Nj).

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

-вычисления значений функции давления Pi(i=Nj-1, Nj-2..0).

-проверка условий итерационного процесса по формуле ( 11). Если итерационный процесс выполняется, то вычисление осуществляется по второму этапу, иначе вышеперечисленные пункты алгоритма вычисления повторяются.

Второй этап алгоритма осуществляется анологичным образом для k+1 временного слоя по направлению переменной y при фиксированном значении x. При этом используется формула (13). Здесь, если итерационный процесс выполняется, то далее вычисление осуществляется по следующему временному слою. Иначе, вышеперечисленные пункты второго этапа повторяются.

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

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

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

Рис. 1. Пользовательский интерфейс программы

Рис. 2.Падение давления наскважине при 1-=0.01, 2-=0.02, 3-=0.03.

Рис. 3. Профили распределения давления при 1-й, 2-й и 3-й годах разработки (=0.02)

Рис. 4. Профили распределения давления при значении:1-=0.01,2-=0.02,3-=0.03.

Вычислительные эксперименты проведены при следующих значениях параматров:длина пласта L=7 км; мощность пласта h=20 м; начальное пластовое давлениеpн=200 атм.; пористость m=0.1; проницаемость пласта 0.01 сп; дебиты скажены 40000 м3/сут. При этом рассмотрены варианты, когда коэффициент вязкости газа принимал значеня =0.01 сп, =0.02 и =0.03. на рис.2 приведены интенсивность падения давления на скважине при различных значениях вязкости газа (=0.01 сп, =0.02 и =0.03) 3-й год разработки. Из рисунка видно, что при большом значении вязкости газападение давления в скважине увеличивается. Распределение давления в пласте приводится в сечении y=0.5 при различных годах разработки (рис. 3) и при различных вариантах значение вязкости газа (рис. 4). Из рис.4 видно показывает, что увеличение значения вязкоси газа приводит к медленному распределению давления в пласте.

Литература:

1. Азиз Х., Сеттари Э. Математическое моделирование пластовых системМосква-Ижевск: Институт компьютерных исследований, 2004. – 416с.
2. Самарский А.А. Введение в теорию разностных схем. Наука. М., 1971.
3. Бельман Р., Калаба Р. Квазилиниаризация и нелинейные краевые задачи. Мир, М., 1968.