Метод ускорения численного решения дифференциального уравнения пьезопроводности модели пласта с двойной пористостью Текст научной статьи по специальности «Математика»

УДК 532.546

БОТ: 10.18698/2309-3684-2024-3-317

Метод ускорения численного решения дифференциального уравнения пьезопроводности модели пласта с двойной пористостью

© Д.Н. Майков1'2, С.С. Макаров2

1ООО "Сиам Мастер", Томская область, Томск, 634003, Россия 2Удмуртский федеральный исследовательский центр Уральского отделения Российской академии наук, Удмуртская Республика, Ижевск, 426067, Россия

В работе представлен метод ускорения численного решения дифференциального уравнения пьезопроводности, на примере описания фильтрации в трещиновато-поро-вом пласте, основанной на модели Уоррена-Рута. Исходная система дифференциальных уравнений, описывающая модель фильтрации от матрицы к трещине, записана через комплексные параметры удельный коэффициент проводимости, долю тре-щинно-кавернозной емкости, и объемную среднюю проницаемость трещин. Предлагаемый метод ускорения численного решения системы дифференциальных уравнений, описывающих модель пласта с двойной пористостью, основан на преобразовании традиционной записи конечно-разностной аппроксимации системы для двух дифференциальных уравнений в одно уравнение. Для получения конечно-разностной аппроксимации параметров использована устойчивая неявная разностная схема. Рассмотрены граничные условия первого и второго рода: граница постоянного давления и непроницаемая граница. Результаты тестовых расчетов по предлагаемому методу сопоставлены с аналитическим решением. При сопоставлении сравнивалось изменение давления в скважине, рассчитанное по численному и аналитическому методу. Давление в скважине рассчитывалось по методу Писмена с определением эффективного радиуса для ячейки сетки Вороного. Проведен численный анализ параметров модели многозабойной скважины в пласте с двойной пористостью с использованием псевдостационарной модели потока. В качестве расчетной сетки использовалась двухмерная декартовая неструктурированная нерегулярная сетка Вороного. Численные расчеты матричных уравнений осуществлялись тремя разными методами: стабилизированный метод бисопряжённых градиентов с ^и(0) предобуславливанием, метод Гаусса-Зейделя с релаксацией, метод Ньютона. Показано, что реализация системы дифференциальных уравнений по предлагаемому методу существенно снижает сложность численного решения и сокращает время расчета моделирования процессов фильтрации и интерпретации параметров при гидродинамическом исследовании скважин.

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

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

которых вводятся упрощения. Принимается, что пласт однородный, процесс фильтрации изотермический, идеальная геометрия скважины, дебит скважины постоянный и т.д. Переход от моделирования процесса добычи с постоянным дебитом к моделированию с изменяющимися дебитами осуществляется с помощью принципа суперпозиции (наложения) по времени [1, 2, 3]. Согласно принципу суперпозиции, сложность расчета давления для каждого последующего дебита растет линейно, а расчет функции всех точек давления целиком имеет квадратичную сложность. Благодаря аналитическим решениям сложность вычислений кратно сокращается, но не всегда полученные решения могут быть пригодны для инженерной практики, так как время расчета сложной модели "скважина - пласт - граница" на сотнях и тысячах режимов работы скважины значительно увеличивается, особенно при решении обратной коэффициентной задачи [4]. В связи с этим, актуальным является построение математических моделей и разработка методов их численного решения.

Так, в работе [5] рассмотрены вопросы математического моделирования трёхмерного ламинарного и турбулентного движения вязкой несжимаемой жидкости в многослойных проницаемых структурах -пористых сетчатых материалах. В работе [6] приведены результаты математического моделирования ламинарной и турбулентной фильтрации жидкой несжимаемой среды в пористых сетчатых материалах. Обзор моделей трещиновато-порового пласта приводится в работе [7]. В работах [8, 9] расчет по модели Уоррена-Рута основан на разностной схеме, состоящей из двух уравнений, описывающих давление в ячейке матрицы и в ячейке трещины. В [8] приводится численное решение модели Уоррена-Рута при помощи программного обеспечения COMSOL. В [9] показан анализ вычислительных затрат по времени при расчете по модели Уоррена-Рута распределения давления в трещиноватом пласте. Оригинальная конечно-разностная схема типа «классики» для модели двойной пористости, позволяющая исключить неизвестное давление в матрице на новом шаге по времени, что дает возможность свести систему уравнений в одно при проведении численного расчета, приведена в [10]. Численное моделирование задач двухфазной фильтрации в трещиновато-пористых средах с использованием модели двойной пористости с сильно неоднородным коэффициентом проницаемости, рассматривается в [11]. В работах [12 - 14] получены оригинальные результаты численного моделирования процесса фильтрации в трещиновато-пористой среде с применением модели двойной пористости.

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

Постановка задачи. За основу примем модельное описание фильтрации в трещиновато-поровом пласте Уоррена-Рута [15]. В модели Уоррена-Рута двойной пористости используется псевдостационарная модель потока PSS (PSS — pseudo steady state), пласт представлен в виде идеализированных геометрических фигур [15, 16]. Система уравнений, описывающая модель фильтрации от матрицы к трещине, по-лученое на основе уравнения пьезопроводности, имеет следующий вид [17]:

f V2pf = ф^ -аkmb(pm -Pf) + q ц dt ц

0 = ФтЬст % + а ^ (pm - pf) dt ц

Здесь kft — объемная средняя проницаемость системы трещин, kmb — объемная средняя проницаемость матрицы, pm — давление в матрице, pf — давление в трещине, с^- — сжимаемость трещины, си — сжимаемость матрицы, ц — вязкость флюида, t — время, а — геометрический коэффициент, 9fb — объемная средняя пористость системы трещин, фтЬ — объемная средняя пористость матрицы, q — иточ-

ник/сток в пластовых условиях.

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

Разработка метода, ускоряющего численное решение. Преобразуем систему уравнений (1), записав её через комплексные параметры удельный коэффициент проводимости Я, долю трещинно-каверноз-ной емкости ю, и объемную среднюю проницаемость системы трещин kft. Удельный коэффициент проводимости Я определяется следующим образом:

ak , r2 4n (n + 2)

Я= mb w , где a = —Ц—, (2)

kfb hm

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

Выразим akmb из уравнения (2):

akmb = ^, (3)

r

w

Общая упругоемкость пласта (ф—) определяется из следующего выражения

(ф—, )ш +' = ФшЪ—ш + Фй—' , (4)

здесь — — общая сжимаемость, т — матрица, / — система трещин.

Доля трещинно-кавернозной емкости с определяется следующим образом:

с= —' . (5)

Фй—' + ФшЪ—ш

Выразим фйс и ФшЬсш из уравнения (5), использовав (4)

Фй,—' = С(ф—, )ш . (6)

Сложим к числителю уравнения (5) и вычтем фшЬсш :

ю = Фйъ—' +ФшЪ —ш ФшЪ —ш ^

Фтсг +ФшЪсш .

Далее разделим дробь на две части и сократим первую часть, получим:

с =!--фшЪ—ш-. (8)

Фй—' +ФшЪ—ш

Из получившегося уравнения выведем фшЬсш :

фшЪ—ш =(1 -с)(ф—, )ш+' . (9)

Перепишем систему уравнений (1) использовав (3), (6) и (9). Имеем запись вида:

7 ^ й =и(ф—)»+' ^" 7 к?(- й)+*

" (10)

0 = (1 -о)(ф—, )ш+' ^ + | '(Рш - Р' ).

Полученная система уравнений (10), записанная через комплексные параметры Я, со, и к^, является предпочтительней при моделировании процессов фильтрации и интерпретации результатов ГДИС, так как не требует задания отдельных параметров пористой матрицы и трещины исследуемой системы.

Запишем конечно-разностную аппроксимацию системы уравнений (10) методом контрольных объемов. Приведем q источник/сток к

поверхностным условиям, умножив и разделив на объемный коэффициент В. Аппроксимацию коэффициента проводимости Т произведем по схеме «вверх по потоку». Для аппроксимации параметров по времени будем использовать устойчивую неявную разностную схему. Конечно-разностная аппроксимация для одномерного случая примет вид:

С (Ры - р )п - г ( : - Р )п=« V; (: -:-1)' -V Г2 ^ (РП - РП)+&'

(11)

О = (1 -ю)(фс1)

РП - Р

|П-1

" Ш1 Ш1

ш+:

Т =

Т 1+0.5

м

: (рт - Р )П,

г. и

(кь Л).

±0.5

Дх

±1

ив

±0.5

Гы =

V (фс )

ш+:

В

(12)

(13)

Здесь Т — коэффициент проводимости на грани, п — индекс времени, 0 — объемный расход на поверхности (дебит), Л — площадь

грани между узловыми точками, V — объем ячейки, кг — относительная проницаемость (при однофазной фильтрации кг = 1), — объемная средняя проницаемость системы трещин по рассматриваемой оси, Дхш — расстояние между расчетными узлами, индекс ±0.5

обозначает грани между расчетными узлами.

Преобразуем уравнение (11). Для этого выразим из второго уравнения Р" .

(1 -Ю)(Фс )

РП = ■

ДI

\ к

ш+: рп-1 , к^рп

1 Ш1 + 2 1 й

г и

(1 -ю)(фС1)

ДI

ш+: + А ^и г2 и

(11)

Подставим уравнение (14) в первое уравнение системы (11), полу-

чим:

Т! (1+1 - Р)" - Т1 (р„-. - Р )п =® т: (рп" - рГ1)

- Г э

в

УРШ-1 +РРП

У + Р

- рп 1 г.

Дг

+0П,

(12)

где

ß=4 f. г-

(1 -ю)(фС()

m+f

At

Уравнение (15) запишем следующим образом:

j-1

Et; (j p )n-« Vj (f - f—1)-V ß

V

/ mi

г+ß

ßPfn - P

pi

Л

(13)

Q (14)

/

где г — рассматриваемая расчетная ячейка, ] — соседняя расчетная ячейка.

Перенесем все неизвестные в уравнении (17) в левую часть, тогда окончательно получим запись

' Ni V

j -1

д

i—

ß

i V

V

г+ß

+ ю-

Vl

At

N

P n _ Tn p n —

p E Tij pfj -

j -1

-ю ^ pn— + V- ß

At

д

Vn—1Л

I mi

г+ß

(15)

—on.

Таким образом, после преобразования получаем систему из N уравнений, а не 2N, где N — количество расчетных узлов.

Начальные и граничные условия. При t - 0, давление во всех ячейках равно начальному давлению p . Граничные условия формируются в соответствии с взаимодействием пласта с его законтурной областью. Рассмотрим граничные условия первого (условие Дирихле) и второго рода (условие Неймана).

Наиболее часто используемым граничным условием первого рода является «граница постоянного давления» (19), второго рода — «непроницаемая граница» (20).

Р(0, t) - const. (16)

сХ

- 0.

(17)

/ x-0

Для того чтобы учесть граничные условия применим метод отражения, суть которого заключается во введении фиктивного расчетного узла. Рассмотрим ячейку под индексом г = 1, которая граничит с фиктивным расчетным узлом под индексом г = 0 . Конечно-разностная аппроксимация для ячейки г = 1 при границе постоянного давления (21) и непроницаемой границе (22) будет иметь следующий вид:

TL (Р2 — pf1 )П — Г (f — pf1 )П - D .

(18)

С (P2 — P )П - D .

(19)

л = ш Ъ (Р«П - Р«П -)-VР

УРШП-1 +РРП У + Р

- рп ± г.

+ &

(20)

Численный анализ применения метода ускорения. Рациональность применения предлагаемого метода ускорения рассмотрим на примере численного решения дифференциального уравнения пьезо-проводности модели пласта с двойной пористостью.

Для численных расчетов построена двухмерная декартовая неструктурированная нерегулярная сетка Вороного. Построение сетки осуществлялось при помощи алгоритма Форчуна [18], так как данный алгоритм обладает сложностью O(nlogn). Локальное измельчение ячеек вокруг скважины осуществлялось согласно положениям, изложенным в работе [19]. Моделирование бесконечного пласта, осуществилось при помощи построения квадратной сетки с длиной стороны м. Время работы скважины составляло 105 ч, при расчетах это время было разбито на 101 расчетный узел в диапазоне с 10-5 до 105 ч. по 10 значений на логарифмический шаг. Для проверки корректности расчетов давления по модели [15] проводилось сопоставление с аналитическим решением модели многозабойной скважины, вскрывающей пласт с двойной пористостью вертикально по всей толщине [20]. Моделирование многозабойной скважины с полным вертикальным вскрытием пласта при численном решении осуществлялось согласно условию, что изменение давления одинаково в каждом стволе скважины [20, 21]. Для этого давление во всех п ячеек, в которых расположены стволы скважины, рассматривались как единая ячейка с притоком, равным общему дебиту скважины. Переток между стволом и расчетной ячейкой, в которой расположен ствол, рассчитывался по методу Писмена [22], где расчет эквивалентного радиуса ячейки Г0 для сетки Вороного определялся согласно положениям, изложенным в работе [23]:

г0 = ехР

!П (4 К

(21)

где I — индекс ячейки скважины, у — индекс соседней ячейки, Ь — длина грани между ячейками I и у , ^ — дистанция между центрами ячеек I и у, # — угол потока.

Матричные уравнения, составленные на основе уравнений (18), (20), решались тремя методами:

• Стабилизированный метод бисопряжённых градиентов (BiCGStab) с ГШ(0) предобуславливанием [24];

• Метод Гаусса-Зейделя с релаксацией [25] (коэффициент релаксации равен 0.91);

• Метод Ньютона [26].

Метод Ньютона использовался без оптимизации разреженных матриц, матрица производных рассчитывалась численно. Для сравнения сложности вычислений конечно-разностная аппроксимация модели двойной пористости PSS решалась методом Ньютона согласно системе уравнений (11), в котором вычислялись два неизвестных параметра: давление в блоке матрицы Рт, давление в блоке трещины

Р . Параметры остановок итерационных алгоритмов подбирались таким образом, чтобы среднее отклонение относительно аналитического решения [12] по рассчитываемому давлению в скважине считается одинаковым для всех и составляло 0.82 %.

Исходные данные для расчетов приведены в таблицах 1 и 2, в качестве граничных условий использовалась «непроницаемая граница».

Таблица 1

Основные параметры модели

Параметры Значения

Толщина пласта, h 10 м

Проницаемость пласта, к 100 мД

Общая сжимаемость, с 0.00005 атм-1

Пористость, <р 0.2

Объемный коэффициент, В 1 м3/ст.м3

Вязкость флюида, л 5 сПз

Общий дебит, q 100 м3/сут

Удельная проводимость, со 0.1

Емкостной коэффициент, X 0.00001

Таблица 2

Параметры ответвлений

Ствол № X, м Y, м Радиус, м

1 0 0 0.065

2 50 50 0.09

3 500 250 0.045

Расчетная сетка показана на рисунке 1. Количество расчетных ячеек составило 1400. Результаты расчетов показаны в таблице 3 на диагностическом графике Бурде [27], приведенном на рисунке 2.

Из таблицы 3 видно, что наиболее быстрым методом решения уравнения (18), (20) является BiCGStab с ГШ{0) предобуславлива-нием, чуть больше времени потребовалось для расчёта по методу Гаусса-Зейделя с релаксацией. Расчет по не оптимизированному методу Ньютона оказался в 50 раз медленнее, чем BiCGStab с ГШ(0)

предобуславливанием и в 46 раз медленнее, чем метод Гаусса-Зейделя

с релаксацией. 52000

41600

31200 20800

10400

0

-10400

-20800 -31200 -41600

-52000

-53638.273 - 42911 - 32184 - 21457 -10730 - 3 10724 21451 32178 42905 53638

Длина, м

Рис. 1. Расчетная сетка

Таблица 3

Результаты расчетов

Метод расчета Решаемые уравнения Время расчетов, с

BiCGStab с ШЩ0) предобуславлива- (18), (20) 1.21

нием

Метод Гаусса-Зейделя с релаксацией (18), (20) 1.33

Метод Ньютона (18), (20) 61.33

Метод Ньютона (11) 121.89

Из расчетов видно, что решение уравнения (11) по методу Ньютона в 2 раза больше по времени, чем решение уравнений (18), (20) записанных по предлагаемому методу ускорения.

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

многозабойной скважины в пласте с двойной пористостью PSS, и аналитическому решению [20].

5

Е 9

- 0 Изменение давления, численная модель „'100 ° Производная изменения давления, численная модель

§ , -Изменение давления, аналитическая модель

д Э — Производная изменения давления, аналитическая модель 2 2

1 10 ..............

5 5

Ь 2

х ^

1 2 | 0,1

I 5

I 10"5 Ю-4 0,001 0,01 0,1 1 10 100 1000 ю4 ю5 £ Время, час

Рис. 2. Изменение давления в многозабойной скважине в пласте с двойной пористостью Г88

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

Для расчета модели с граничными условиями воспользуемся данными таблицы 1 и таблицы 2, а также построим новую расчетную сетку со стороной 1500 м. Расчеты производились при двух границах: все четыре стороны — непроницаемая граница и все четыре стороны — граница постоянного давления. Результаты расчетов показаны на рисунке 3.

1000

3

¡100

Изменение давления, непроницаемая граница Производная изменения давления, непроницаемая граница — Изменение давления, граница постоянного давления ^ Производная изменения давления, граница постоянного давления

....................................................оо0

10

5 х

2

¡0,1

10~5 Ю^4 0,001 0,01 0,1 1 10 100 1000 104 105

Время, час

5

5

2

5

Рис. 3. Изменение давления в многозабойной скважине в пласте с двойной пористостью PSS при заданных граничных условиях

На рисунке 3 видно, что возмущение от скважины достигает границ за 600 часов, после чего в случае с границей постоянного давления перепад становится постоянным, а при непроницаемой границе — резко возрастает.

Выводы. Предложен новый метод ускорения численного решения уравнения пьезопроводности модели пласта с двойной пористостью. Численное решение по предлагаемому методу согласуется с аналитическим решением. Запись конечно-разностной аппроксимации системы дифференциальных уравнений по предлагаемому методу, по сравнению с традиционным методом, снижает сложность численного решения и существенно сокращает время при сохранении точности расчета. Применение предлагаемого метода ускорения расчета дифференциального уравнения пьезопроводности позволит оперативно производить интерпретацию ГДИС с учетом комплексных параметров трещиновато-порового коллектора системы "скважина-пласт".

ЛИТЕРАТУРА

[1] Walker A.C. Estimating reservoir pressure using the principle of superposition. SPE, 1968. DOI: 10.2118/2324-MS

[2] Cinco-Ley H., Samaniego F. Use and misuse of the superposition time functionin well test analysis. SPE-19817,1989. DOI: 10.2118/19817-MS

[3] Майков Д.Н., Исупов С.В., Макаров С.С., Аниканов А.С. Метод ускорения расчета давления при изменяющихся дебитах по истории эксплуатации скважины. Нефтяное хозяйство, № 9, 2021, с. 105-107. DOI: 10.24887/00282448-2021-9-105-107

[4] Майков Д.Н., Макаров С.С. Численное исследование алгоритмов оптимизации при адаптации гидродинамической модели по результатам исследований скважин. Математическое моделирование, 2022, том 34, № 9, с. 71-82.

[5] Городнов А.О., Лаптев И.В., Сидоренко Н.Ю. Математическое моделирование процессов ламинарной и турбулентной фильтрации жидкой несжимаемой среды в пористых сетчатых материалах. Математическое моделирование и численные методы, 2023, № 2, с. 67-89. DOI: 10.18698/2309-3684-20232-6789

[6] Холмуродов А.Э., Дильмурадов Н. Математическое моделирование одномерного нелинейного движения в насыщенной жидкостью пористой среде. Математическое моделирование и численные методы, 2018, № 1, с. 3-15. DOI: 10.18698/2309-3684-2018-1-315

[7] Блонский А.В., Митрушкин Д.А., Савенков Е.Б. Моделирование течений в дискретной системе трещин: физико-математическая модель. Препринты ИПМим. М.В. Келдыша, 2017, № 65, с. 28. DOI: 10.20948/prepr-2017-65

[8] Maaref S., Moosavi R., Riazi M. Comparison of numerical methods to solve Warren-Root PDE of naturally fractured reservoirs. The 1st National Conference on Oil and Gas Fields Development, 2015.

[9] Истрафилов М.Я., Морозкин Н.Н. Решение задачи фильтрации в трещиноватом пласте с использованием модели Уоррена-Рута. Вестник Башкирского университета, 2017, Т. 22, № 1, с.15-19.

[10] Афанаскин И.В., Вольпин С.Г., Ломакина О.В., Штейнберг Ю.М. Интерпретация межскважинных исследований карбонатных коллекторов методом двух режимов с помощью численных моделей. Программные продукты и системы, 2018, № 3, с. 500-506.

[11] Васильев В.И., Васильева М.В., Григорьев А.В., Прокопьев Г.А. Математическое моделирование задачи двухфазной фильтрации в неоднородных трещиновато-пористых средах с использованием модели двойной пористости и метода конечных элементов. Ученые записки Казанского университета. Серия: Физико-математические науки, 2018, Т. 160, № 1, с. 165-182

[12] Вабищевич П.Н., Григорьев А.В. Численное моделирование фильтрации флюида в анизотропной трещиновато-пористой среде. Сибирский журнал вычислительной математики, 2016, Т. 19, № 1, с. 61-74. DOI: 10.15372/SJNM20160105.

[13] Вабищевич П.Н., ГригорьевА.В. Численное моделирование фильтрации на основе модели двойной пористости. Суперкомпьютерные технологии математического моделирования: труды второй международной конференции, 2013, с. 50-61.

[14] Григорьев А.В. Численное моделирование фильтрации в трещиновато-пористой среде. Математические заметки ЯГУ, 2013, Т. 20, № 2, с. 237-245.

[15] Warren J.E., Root P.J. The behavior of naturally fractured reservoirs. SPE Journal, 1963, pp. 244-255. DOI: 10.2118/426-PA

[16] Stewart G. Well test design and analysis. PennWell Corp., 2011, pp. 549-603.

[17] Abdus Satter, Ghulam M. Iqbal. Reservoir Engineering: The Fundamentals, Simulation, and Management of Conventional and Unconventional Recoveries. Texas, Gulf Professional Publishing, 2016, 486 p.

[18] Fortune S. A sweepline algorithm for Voronoi Diagrams. Algorithmica, 1987, vol.

2, pp. 153-174.

[19] Zhang L., Wenshu Zha, Lu D. Generation and application of adaptive PEBI grid for numerical well testing (NWT). Proceedings 2013 International Conference on Mechatronic Sciences, 2013. DOI: http://dx.doi.org/10.1109/MEC.2013.6885542

[20] Майков Д.Н., Макаров С.С. Параметрический анализ модели многозабойной скважины в пласте с двойной пористостью. Вестник ТюмГУ. Физико-математическое моделирование. Нефть, газ, энергетика, 2023, Том 9, № 3 (35), с. 100-116.

[21] Майков Д.Н., Борхович С.Ю. Аналитическая модель многозабойной скважины с полным вертикальным вскрытием пласта. Нефть. Газ. Новации, 2020, № 11 (240), с. 61-65.

[22] Peaceman D.W. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. Society of Petroleum Engineers Journal, 1983, 23(03), pp. 531-543.

[23] Syihab Z. Simulation on Discrete Fracture Network Using Flexible Voronoi Grid System. Petroleum Engineering, Texas A&M University, 2009, 194 p.

[24] Saad Y. Iterative methods for sparse linear systems. SIAM, 2003, 184 p.

[25] Глазырина Л.Л., Карчевский М.М. Введение в численные методы: учебное пособие. Казань, Казанский университет, 2017, с. 122.

[26] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Москва, БИНОМ. Лаб. знаний, 2011, 636 с.

[27] Bourdet D. et al. A new set of type curves simplifies well test analysis. World Oil, 1983, pp. 95-106.

Статья поступила в редакцию 01.04.2024

Ссылку на эту статью просим оформлять следующим образом:

Майков Д.Н., Макаров С.С. Метод ускорения численного решения дифференциального уравнения пьезопроводности модели пласта с двойной пористостью. Математическое моделирование и численные методы, 2024, № 3, с. 3-17.

Майков Дмитрий Николаевич — ведущий специалист, ООО "СиамМастер", младший научный сотрудник ФГБУН УдмФИЦ УрО РАН. e-mail: [email protected]

Макаров Сергей Сергеевич — д-р техн. наук, ведущий научный сотрудник, ФГБУН УдмФИЦ УрО РАН. e-mail: [email protected]

The efficient method for numerical solution of dual porosity reservoir model

© D.N. Maykov1'2, S.S. Makarov2

1SIAM MASTER Ltd, Tomsk region, Tomsk, 634003, Russia

2Udmurt Federal Research Center of the Ural Branch of the Russian Academy of Sciences, Udmurt Republic, Izhevsk, 426067, Russia

This paper presents a method for accelerating the numerical solution of the diffusivity equation based on the Warren-Root model. A general differential equation system describing the filtration model from the matrix to the fracture is written through the complex parameters of the transmissivity ratio, the storativity ratio, and the volumetric average permeability of the fracture system. The proposed method for accelerating the numerical solution of a differential equation system describing a double-porosity reservoir model is based on converting the traditional form of a finite-difference approximation of the system for two differential equations into one equation. A stable implicit difference scheme is used to obtain a finite-difference approximation of the parameters. Boundary conditions of the first and second kind are considered: a constant pressure boundary and an impermeable boundary. The results of the test calculations using the proposed method are compared with the analytical solution. The pressure change in the well was compared, calculated by numerical and analytical methods. The pressure in the well was calculated using the Peaceman method with the effective radius for the Voronoi grid cell. A numerical analysis of the parameters of a multilateral well in a double porosity formation model is carried out. A two-dimensional Cartesian unstructured irregular Voronoi grid was used as the calculated grid. The numerical calculations of the matrix equations were carried out by three different methods: the biconjugate gradient stabilized method with ILU(0) preconditioning, the Gauss-Seidel method with relaxation, and the Newton method. It is shown that the implementation of a differential equation system according to the proposed method significantly reduces the complexity of the numerical solution and reduces the calculation time of the filtration process modeling and pressure transient analysis interpretation.

Keywords: model, dual porosity reservoir numerical solution, multilateral well, diffusivity equation, parametric analysis

REFERENCES

[1] Walker A.C. Estimating reservoir pressure using the principle of superposition. SPE, 1968. DOI: 10.2118/2324-MS

[2] Cinco-Ley H., Samaniego F. Use and misuse of the superposition time functionin well test analysis. SPE-19817,1989. DOI: 10.2118/19817-MS

[3] Maykov D.N., Isupov S.V., Makarov S.S., Anikanov A.S. Metod uskoreniya rascheta davleniya pri izmenyayushchihsya debitah po istorii ekspluatacii skva-zhiny [The efficient method for pressure calculation at variable rate]. Neftyanoe hozyajstvo [Oil Industry], № 9, 2021, pp. 105-107. DOI: 10.24887/0028-24482021-9-105-107

[4] Maykov D.N., Makarov S.S. Numerical investigation of optimization algorithms for the hydrodynamic model adaptation based on the well test results. Mathematical Models and Computer Simulations, 2022, vol. 34, № 9, pp. 71-82.

[5] Gordonov A.O., Laptev I.V., Sidorenko N.Yu. Mathematical modeling of laminar and turbulent filtration processes of liquid in compressible medium in porous meshmaterials. Mathematical modeling and computational methods, 2023, № 2, pp. 67-89. DOI: 10.18698/2309-3684-2023-2-6789

[6] Kholmurodov A.E., Dilmuradov N. Mathematical modeling of one-dimensional non-linear motion in a fluid-saturated porous medium. Mathematical modeling and computational methods, 2018, № 1, pp. 3-15. DOI: 10.18698/2309-36842018-1-315

[7] Blonsky A.V., Mitrushkin D.A., Savenkov E.B. Discrete fracture network modelling: physical and mathematical model. Keldysh Institute PREPRINTS, 2017, № 65, p. 28. DOI 10.20948/prepr-2017-65

[8] Maaref S., Moosavi R., Riazi M. Comparison of numerical methods to solve Warren-Root PDE of naturally fractured reservoirs. The 1st National Conference on Oil and Gas Fields Development, 2015.

[9] Istrafilov M.Y., Morozkin N.N. Solution of the filtration problem in the fissured formation using warren-root model. Vestnik Bashkirskogo universiteta, 2017, vol. 22, no. 1, pp. 15-19.

[10] Afanaskin I.V., Volpin S.G., Lomakina O.V., ShteynbergYu. M. Carbonate reservoirs crosswell survey interpretation by a two-rate test using numerical models. Software & Systems, 2018, no. 3, pp. 500-506.

[11] Vasilyev V.I., Vasilyeva M.V., Grigorev A.V., Prokopiev G.A. Mathematical modeling of the two-phase fluid flow in inhomogeneous fractured porous media using the double porosity model and finite element method. Lobachevskii Journal of Mathematics, 2018, vol. 160, no. 1, pp. 165-182.

[12] Vabishchevich P.N., Grigoriev A.V. Numerical modeling of a fluid flow in anisotropic fractured porous media. Numerical Analysis and Applications, 2016, vol. 19, no. 1, pp. 61-74. DOI: 10.15372/SJNM20160105.

[13] Vabishchevich P.N., Grigoriev A.V. CHislennoe modelirovanie fil'tracii na osnove modeli dvojnoj poristosti [Numerical simulation of filtration based on the dual porosity model]. Superkomp'yuternye tekhnologii matematicheskogo modelirovaniya: trudy vtoroj mezhdunarodnoj konferencii [Second international conference Supercomputer technologies of mathematical modelling], 2013, pp. 50-61.

[14] Grigoriev A.V. Chislennoe modelirovanie filtracii v treshchinovato-poristoj srede [Numerical modeling of filtration in a fractured porous medium].

Matematicheskie zametki YAGU [Mathematical Notes of NEFU], 2013, vol. 20, no. 2, pp. 237-245.

[15] Warren J.E., Root P.J. The behavior of naturally fractured reservoirs. SPE Journal, 1963, pp. 244-255. DOI: 10.2118/426-PA

[16] Stewart G. Well test design and analysis. PennWell Corp., 2011, pp. 549-603.

[17] Abdus Satter, Ghulam M. Iqbal. Reservoir Engineering: The Fundamentals, Simulation, and Management of Conventional and Unconventional Recoveries. Texas, Gulf Professional Publishing, 2016, 486 p.

[18] Fortune S. A sweepline algorithm for Voronoi Diagrams. Algorithmica, 1987, vol. 2, pp. 153-174.

[19] Zhang L., Wenshu Zha, Lu D. Generation and application of adaptive PEBI grid for numerical well testing (NWT). Proceedings 2013 International Conference on Mechatronic Sciences, 2013. DOI: http://dx.doi.org/10.1109/MEC.2013.6885542

[20] Maykov D.N., Makarov S.S. Parametricheskij analiz modeli mnogozabojnoj skvazhiny v plaste s dvojnoj poristost'yu [The multilateral well model with complete vertical opening of naturally fractured reservoir]. Vestnik TyumGU. Fiziko-matematicheskoe modelirovanie. Neft', gaz, energetika [Tyumen State University Herald. Physical and Mathematical Modeling. Oil, Gas, Energy], 2023, vol. 9, no. 3(35), pp. 100-116.

[21] Maykov D.N., Borkhovich S.Yu. Analiticheskaya model' mnogozabojnoj skvazhiny s polnym vertikal'nym vskrytiem plasta [Analytical model of a multilateral well with full vertical penetration of the formation]. Neft'. Gaz. Novacii [Oil. Gas. Innovations], no. 11 (240), 2020, pp. 61-65.

[22] Peaceman D.W. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. Society of Petroleum Engineers Journal, 1983, 23(03), pp. 531-543.

[23] Syihab Z. Simulation on Discrete Fracture Network Using Flexible Voronoi Grid System. Petroleum Engineering, Texas A&M University, 2009, 194 p.

[24] Saad Y. Iterative methods for sparse linear systems. SIAM, 2003, 184 p.

[25] Glazyrina L.L., Karchevskij M.M. Vvedenie v chislennye metody: uchebnoe posobie [Introduction to numerical methods: textbook]. Kazan', Kazanskij uni-versitet, 2017, 122 p.

[26] Bahvalov N.S., Zhidkov N.P., Kobelkov G.M. CHislennye metody [Numerical methods]. Moscow: BINOM. Lab. znanij [BINOM. Lab. knowledge], 2007, 636 c.

[27] Bourdet D. et al. A new set of type curves simplifies well test analysis. World Oil, 1983, pp. 95-106.

Maykov D.N., Leading Specialist SIAM MASTER Ltd, Junior Researcher Udmurt Federal Research Center of the Ural Branch of the Russian Academy of Sciences. e-mail: [email protected]

Makarov S.S., Dr. Sci. (Eng.), Leading Researcher, Udmurt Federal Research Center of the Ural Branch of the Russian Academy of Sciences. e-mail: [email protected]