______ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА________________
2008 Математика. Механика. Информатика Вып. 4(20)
УДК 532.546
Модель течения жидкости в дискретной сети сжимаемых трещин
А. А. Щипанов, С. В. Русаков
Пермский государственный университет, 614990, Пермь, ул. Букирева, 15
Предложена новая модель дискретной сети сжимаемых трещин для описания течения жидкости в деформируемой горной породе. Рассматривается течение однородной жидкости в сети вертикальных пересекающихся трещин, рассекающих непроницаемую горную породу на блоки, которые составляют матрицу. Связь раскрытости трещины и ее проницаемости введена исходя из сопоставления решения Пуазеля и закона Дарси. Эволюция давления в процессе нестационарного течения жидкости вызывает изменение раскрытости трещин, обусловленное сжимаемостью матрицы (горной породы). Предложенная математическая модель составляет основу численной процедуры моделирования течения в сети сжимаемых трещин. Для двумерного случая на основе интегро-интерполяционного метода построена конечно-разностная схема и разработан пакет программ, позволяющий моделировать течение жидкости в сети сжимаемых трещин к источнику или стоку (добывающей или нагнетательной скважине). Разработанная модель и пакет программ имеют практическое приложение при оценке параметров деформируемых трещиноватых пластов месторождений углеводородов, в частности, при интерпретации гидродинамических исследований скважин.
Введение
К настоящему времени разработан комплекс математических моделей и технологий, позволяющих моделировать течения в природных пластах, подверженных естественной трещиноватости, с учетом всего многообразия механизмов извлечения углеводородов [1, 2]. Одной из основных проблем практического применения разработанных технологий является идентификация параметров систем трещин [3-6]. Возможным путем решения данной проблемы является разработка теоретических подходов, основанных на использовании аналитических решений и численных методов, которые позволяют выполнить идентификацию свойств трещин посредством интерпретации эмпирических данных. Практически
© А. А. Щипанов, С. В. Русаков, 2008 Исследования выполнены при поддержке CRDF, грант Y3-MP-09-07 и РФФИ, грант "р_офи" №0701-97618
важной представляется разработка подходов, позволяющих оценить весь комплекс свойств системы трещин.
Наиболее широкое распространение для решения проблемы идентификации трещиноватости в природных пластах в последнее время получила концепция дискретной сети трещин (DFN - Discrete Fracture Network). На основе данного подхода производится описание сложной геометрии природной трещиноватости и моделирование течения в системе трещин. Предложенная для описания течения в трещиноватых массивах горной породы в приложении к задачам гидрогеологии [3], концепция была расширена для решения проблем разработки месторождений углеводородов [4-6]. В работах [5, 6] отражено ее развитие для моделирования системной трещиноватости, коридоров трещин и сейсмических нарушений. Моделирование течения в дискретной сети трещин позволяет решать задачи теории перколяции [7]. Основное внимание в
этой теории, которая находит практическое приложение при идентификации и моделировании природных пластов, уделяется связанности сети трещин.
Разработанная на основе исследований [3-6] методология позволяет описать реальную геометрию трещин, выделить семейства трещин, выполнить построение синтетической системы трещин и осуществить моделирование течения в этой системе, что, например, позволяет произвести оценку параметров системы трещин по данным промысловых исследований. В результате практического применения разработанной методологии возможно построение модели системы трещин на уровне пласта и вычисление их эквивалентных параметров на дискретной пространственной сетке, которая используется для моделирования многофазных потоков, происходящих в пластах месторождений углеводородов при их разработке.
Рассмотренные выше исследования предполагают учет влияния сжимаемости трещин только в зависимости пористости от давления, при этом проницаемость системы трещин предполагается постоянной. Анализ эмпирических данных показывает, что проницаемость трещин в более существенной степени, чем пористость, зависит от изменения пластового давления [8-9]. Учет этого аспекта необходим при моделировании течений в природной системе трещин.
Идентификация атрибутов трещин (плотность, раскрытость, ориентация, падение) и континуальных свойств трещиноватого пласта (пористость, проницаемость системы трещин, размер блоков матрицы, на которые разбивается горный массив этой системой) обычно выполняется путем интерпретации эмпирических данных (сканирование стенки скважин, исследование керна, гидродинамические исследования скважин, история разработки залежи и др.). При этом введение однозначных взаимосвязей между атрибутами трещин и свойствами системы трещин (например, между раскрытостью трещин и пористостью, проницаемостью) в рамках предложенных подходов не предполагается, что обусловлено трудностями в оценке атрибутов трещин. Интерпретация результатов современных методов исследований пласта позволяет оценить атрибуты трещин, определяю-
щие поток (длину, связность, проводимость и, как результат, проницаемость трещин). Оценка емкостных свойств (раскрытость, пористость) вызывает затруднения.
Идеализированный подход, когда отдельная трещина представляется пустотным пространством между двумя параллельными пластинами (границами горной породы) позволяет связать пористость и проницаемость с раскрытостью трещины на основе сопоставления закона Дарси и решения для течения Пуазеля и широко используется при решении прикладных задач разработки природных пластов [10]. Допущение линейности стенок трещин является существенной идеализацией, поскольку в природных пластах они обладают высокой шероховатостью - в связи с этим часто переходят к понятию гидравлической апертуры (раскрытости). Однако безусловным преимуществом данного подхода является связь фильтрационных (проницаемость) и емкостных (пористость) свойств трещин с ее раскрытостью, что позволяет успешно решать задачи идентификации всего комплекса свойств трещин.
Представленный обзор делает актуальной и практически важной задачу разработки теоретического подхода, описывающего сеть трещин как систему, в которой все атрибуты элементов и свойства сети в целом имеют определенную взаимосвязь. Основу разрабатываемого подхода может составлять современная и широко применяемая на практике концепция дискретной сети трещин (DFN). Исходя из анализа эмпирических данных, необходимо развитие данной концепции для комплексного учета сжимаемости горной породы (матрицы) и трещин и ее влияния на континуальные свойства системы трещин.
1. Модель дискретной сети
сжимаемых трещин
Рассматривается течение однородной жидкости в сети вертикальных пересекающихся трещин, рассекающих горную породу (предполагается, что порода непроницаемая для жидкости) на блоки, которые составляют матрицу. Сеть трещин состоит из двух семейств с одинаковыми значениями атрибутов (плотность, раскрытость) и ортогональными направлениями. Двумерное распределение дискретных элементов позволяет описать
представленную сеть сжимаемых трещин (CDFN - Compressible Discrete Fracture Network). Основными атрибутами элемента сети являются линейный размер1 L и раскры-тость трещин b . Дискретный элемент сети сжимаемых трещин приведен на рис. 1.
Линейный размер элемента складывается из линейного размера блока матрицы и раскрытости трещины:
L=L„+b = H:+b{)
(1)
1 dVr >r dp
(2)
vr =i;uexp Kcrb-p
Введем коэффициент изменения линейного размера блока матрицы:
Lr = Z"exp lcL (5)
Сравнивая соотношения (4) и (5), получим:
cL =cr/2.
Представленный подход позволяет выразить зависимость раскрытости трещины от распирающего давления жидкости, используя коэффициент изменения размера блока матрицы (на основе раскрытости и размера блока матрицы при некоторых определенных условиях (опорном давлении)). Из соотношений (1) и (5) получаем:
Рис. 1. Элемент дискретной сети сжимаемых трещин (СОПЧ-элемент)
Сжимаемость матрицы (горной породы) определяется следующим образом:
Интегрируя данное соотношение, можно получить зависимость объема блока матрицы от распирающего давления жидкости в трещине:
(3)
b = b + Lr 1-exp ♦ cL p-p =
= Ь-Ё], exp 4-cL
Пористость элемента сети выражается как ф=§ iL-b ~Jj} .
Проницаемость каждой трещины2, которая входит в элемент, получаем по аналогии с подходом [10] на основе сопоставления решения для течения Пуазеля и закона Дарси:
kx=kv = b3 ^-Ы2~]}21} .
Проницаемость элемента можно оценить как сумму проницаемостей двух пересекающихся трещин:
к = кх +kv = b2<f>j\2 .
Зависимость раскрытости от давления жидкости (6) приводит к зависимости пористости и проницаемости трещин от давления:
В предположении недеформируемости и непроницаемости верхней и нижней границ пласта толщины Н , содержащего рассматриваемую сеть трещин, соотношение (3) можно записать следующим образом:
Фр = г)р tl- bp~~]l2 к р =b2 p~(j> р ~]\2 .
(7)
(8)
V. = н 4., 2 = н t!
exp Ясг
(4)
о
Сг =
2
О
1 Плотность трещин обратно пропорциональна линейному размеру
2 Ориентируем оси координат в соответствии с на-
правлениями трещин
2. Численное моделирование течения в дискретной сети сжимаемых трещин
Рассмотрим двумерную дискретную сеть трещин. Эволюция давления в сети трещин можно описать уравнением3, аналогичным уравнению для пористой среды [11]:
{/иВ ] ді
Ф.
В
+ /•
(9)
а (у-уг
др др\ V
К
V
1 дК
V г
У г Ф
= *-ф£г
исходное уравнение для эволюции давления (9) можно записать в форме
1 -ф ф
с„ +— с
в
В
др
дії
+ /•
(10)
Рассмотрим двумерную сеть трещин. Множество элементов (ячеек разностной сетки) может быть описано координатами их центров (узлов). Таким образом, получаем двумерную сетку:
0,х = х]=1к,у. = уЛ;/,7 = 0,N;Ик=Ьп .
Эволюция давления моделируется на множестве последовательных временных слоев (дис-кретнзацня уравнения по времени):
Преобразуя производную по времени в правой части уравнения (9), учитывая (2), сжимаемость жидкости с [11], а также выражение производной от пористости:
О = =пт,п = 0,М;Мт = Тп
Конечно-разностное уравнение, которое является дискретным аналогом уравнения в частных производных (10) выглядит следующим образом:
тп Хп+1 пп+1 тп Хп+1 пп+1 I
Тх,г+1/2,у Р1+1,у рг,у _ Т х, /—1/2,у Р1,у Р 1 • +
1,3 -
Тп Рп+1 - рп-1 - Тп Р
Ту,і,]-1/2 л?,3 +1 рі,J _ /і7-1/2 л-
п рп+1 _ Рп+1 у,і,7-1/2 л?,3 рі,
43-1
Уравнение (10) в двумерном случае имеет вид
д Г К ф] 5 [ ку ф|
і ______________
дх [/иВ дх\ ду\/иВ ду
1 -ф ф
—-сг+—с В г В
др
дґ
+ /•
(11)
= К
і-ф; ф;
----- с + — с
вп в:
рГ ~рп
У+^+1. (15)
Связанность Т вычисляется с использованием гармонического осреднения проницаемости и осреднения "вверх по потоку" [11] свойств жидкости (произведения вязкости и объемного коэффициента):
Для решения данного уравнения требуется задание начального условия:
грп
Іх,і± 1/2,у
-ІІІ/2,7 ■
р *,у,0=р°;
(12)
используются два типа граничных условий заданное постоянное давление:
Р*, У, *Ь=Рп,
(13)
Коэффициенты нелинейного уравнения (15) являются функциями от давления (в частности (7), (8)), поэтому используются последовательные итерации по нелинейности:
Т п+1,5 П+1,5+1_ ТГ гп+1 .
Р ~
отсутствие потока на границе (непроницаемая граница):
др
дп
= 0.
(14)
Уравнение получено комбинацией законов сохранения массы и Дарси
4 Сжимаемость жидкости определяет зависимость объемного коэффициента (плотности) от давления
+к,
1 -фг
11+1,5
Фг
11 + 1,5
В
———6' + —
«+1,х г
11+1,5
Для решения системы линейных алгебраических уравнений (СЛАУ), полученной исходя из (15), используются итерационные методы Якоби, Гаусса-Зейделя и последовательной верхней релаксации (ПВР) [11].
с
Описанный конечно-разностный метод с двумя типами граничных условий (13), (14) и аналогичными типами условий5 на источнике (стоке), реализован в виде пакета программ, который позволяет моделировать поток и источнику или стоку (скважине, рис. 2) на основе предложенной модели дискретной сети сжимаемых трещин.
3. Влияние сжимаемости трещин на характеристики течения
Новизна предложенного подхода дискретной сети заключается в учете сжимаемости трещин, которая влияет как на пористость, так и на проницаемость. Серия вычислительных экспериментов выполнена для оценки влияния сжимаемости трещин, определяемой сжимаемостью породы (твердой составляющей пласта), на основные характеристики течения жидкости в сети трещин.
Моделируется течение жидкости в квадратной расчетной области к стоку (рис. 2). На границе расчетной области накладывается условие заданного давления (13), равного начальному (12). Сток с постоянной интенсивностью присутствует в первой половине расчетного интервала времени, что вызывает конусообразное понижение давления в расчетной области (наибольшее снижение давления происходит в элементе, содержащем сток, - центральном элементе, рис. 2). Затем сток отключается, что вызывает конусообразное повышение давления в расчетной области. Анализ результатов проведем на базе эволюции давления в элементе, содержащем сток.
Проанализируем результаты расчетов двух вариантов (густота трещин Ь= 1 м-1): сжимаемость породы сг составляет 1.0*10-5 МПа-1 - такую породу можно охарактеризовать как относительно слабосжимаемую; сжимаемость породы составляет 1.0*10-4 МПа-1 - сжимаемая порода. Изменение раскрытости (апертуры) трещин при изменении давления жидкости, а также вызванное этим изменение пористости и проницаемости трещин приведено для второго варианта на рис.3. Аналитическая зависимость построена по выражениям (7)-(8) в заданном интервале давлений, численная - при численном решении уравнения давления (15) и
5 Условия заданного расхода (дебита или приемистости) или заданного давления на стенке скважины (стока или источника)
вычислении зависящих от него параметров (6)-(8).
В случае слабосжимаемой породы понижение и повышение давления происходит монотонно, при этом скорость падения (восстановления) давления уменьшается со временем при включении (отключении) стока (рис. 4).
В случае сжимаемой породы (рис. 5) понижение и повышение давления происходит также монотонно. Однако скорость падения (восстановления) давления сначала снижается, а затем начинает возрастать, что связано с изменением апертуры и соответственно проницаемости трещин. В случае понижения давления снижение проницаемости трещин вызывает необходимость сравнительно большего понижения давления (увеличения депрессии) для обеспечения той же интенсивности стока. Замедление процесса восстановления давления в начальный период времени связано с низкой проницаемостью трещин, обусловленной предшествующим падением давления.
Для анализа процесса восстановления (падения) давления наиболее информативными являются кривые, отражающие динамику производной давления по времени. В практике идентификации природных пластов месторождений углеводородов при интерпретации кривых восстановления (падения) давления, полученных в процессе гидродинамических исследований скважин, используется кривая произведения времени на производную от давления по времени [10]. Анализ кривой позволяет выявить диагностические признаки присутствия непроницаемых границ, высокопроницаемых включений (например, трещин) и других особенностей околоскважинной зоны пористого пласта. Выполненные вычислительные эксперименты позволяют сделать вывод о возможности использования данной кривой для оценки изменения проницаемости трещин.
На рис. 6-9 приведена описанная выше кривая для процесса восстановления давления после отключения стока, на рис. 8-9 использована логарифмическая шкала времени и давления. Немонотонная кривая для слабосжимаемой породы лежит ниже кривой изменения давления (рис. 6, 8). Кривая для сжимаемой породы также немонотонна, но в начальный период времени лежит выше кривой изменения давления (рис. 7, 9). Хорошо видно ускорение восстановления давления, сменяющееся замедлением в этом случае (рис. 7, 9), в то время как для слабосжи-маемой породы характерно замедление (рис. 6, 8). Присутствие ускорения вызвано процессом
восстановления проницаемости сети трещин при восстановлении давления.
Результаты вычислительных экспериментов позволяют сделать вывод о существенном влиянии сжимаемости трещин (определяемой сжимаемостью породы) на характеристики течения в системе трещин. Выявленное ускорение восстановления давления может служить диагностическим признаком присутствия существенного изменения проницаемости среды, который может быть использован при интерпретации гидродинамических исследований скважин.
Ч.80Е+1
Рис. 2. Распределение давления Р в дискретной сети сжимаемых трещин при течении жидкости к стоку (скважине)
time, days
Рис. 4. Эволюция давления (Pressure) в элементе сети (CDFN) после включения (time = 0 дней) и отключения (time = 5 дней) стока; сжимаемость породы составляет 1.0*10-5 МПа'1
Analitical Aperture ■ Numerical Aperture
Pressure, MPa
-------Analitical Porosity --------------Analitical Permeability
■ Numerical Porosity ■ Numerical Permeability
Рис. 3. Зависимость раскрытости (Aperture), пористости (Porosity) и проницаемости (Permeability) трещин от давления (Pressure)
time, days
Рис. 5. Эволюция давления (Pressure) в элементе сети (CDFN) после включения (time = 0 дней) и отключения (time = 5 дней) стока; сжимаемость породы составляет 1.0*10~4 МПа-1
і -1_геІ, сіауз
Рис. 6. Восстановление давления Р после отключения стока. Кривые приращения давления \P-P_ref и произведения времени на производную от давления по времени \(-t_ref\*dP/dt. Сжимаемость породы составляет 1.0*10"5 МПа"1
Рис. 8. Кривые приращения давления \Р-P_ref и произведения времени, на производную от давления по времени ^-t_ref\*dP/dt (рис. 6) в логарифмической шкале времени и давления. Сжимаемость породы составляет 1.0*10-5МПа"1
Рис. 7. Восстановление давления Р после отключения стока. Кривые приращения давления \P-P_ref и произведения времени, на производную от давления по времени \ и t_ref\*dP/dt. Сжимаемость породы составляет 1.0*10-4 МПа"4
Рис. 9. Кривые приращения давления \Р-P_ref и произведения времени, на производную от давления по времени \и t_ref\*dP/dt (рис. 7) в логарифмической шкале времени и давления. Сжимаемость породы составляет 1.0*10-4 МПа-1
Заключение
Концепция дискретной сети трещин развита для комплексного учета сжимаемости - предложен подход дискретной сети сжимаемых трещин, который характеризуется следующим:
a) введена связь между сжимаемостью матрицы (породы), которая определяется с достаточной точностью при исследованиях керна, и сжимаемостью трещин, которая обычно используется при моделировании гидродинамики пласта;
b) вычисление проницаемости и пористости производится через раскрытость трещин, что позволяет связать данные параметры и определить их зависимость от давления;
c) использовано естественное для модели двойной пористости представление Уоррена-Рута [12] матрицы и трещин, что позволяет легко вычислить размеры блоков матрицы и геометрический параметр [1, 2] модели двойной среды.
Вычислительные эксперименты показали существенное влияние сжимаемости трещин на характеристики течения в дискретной сети. Динамика изменения давления отличается существенным образом в случае сжимаемости трещин: присутствует интервал ускорения (замедления) изменения давления на кривой восстановления (падения) давления. Данный аспект имеет практическое значение при интерпретации гидродинамических исследований скважин, дренирующих трещиноватые природные пласты, и может быть использован как диагностический признак присутствия существенно сжимаемой системы трещин.
Предложенный подход дискретной сети сжимаемых трещин может служить основой новых методов идентификации трещиноватых природных пластов. Разработанный пакет программ для моделирования течения в сети трещин позволяет выполнять калибровку атрибутов трещин путем моделирования гидродинамических исследований скважин. Начальное приближение для атрибутов трещин (оценка значений атрибутов) может быть получено, например, с использованием метода идентификации, основанного на предложенном подходе.
Условные обозначения
Ь - раскрытость трещины, м;
В - объемный коэффициент жидкости, без-разм.;
с - сжимаемость жидкости, МПа-1; сг - сжимаемость породы (матрицы), МПа-1; с - коэффициент изменения линейного размера блока матрицы, МПа-1;
/ - источниковый член, с- (сут-);
И - размер ячейки (элемента сети), м;
Н - толщина пласта, м; к - проницаемость трещин, мкм2 (Д), 1 Д ~ 0.987 мкм2;
кх , ку - проницаемость трещин по направлениям, мкм2 (Д);
Ь - линейный размер элемента сети, м;
Ьг - линейный размер матрицы, м;
Ьа - размер области моделирования, м М - количество слоев по времени;
N - количество узлов пространственной разностной сетки (элементов сети); р - давление жидкости, МПа;
ра - давление на внешней границе, МПа;
£ - время, с (сут);
Т - связанность трещин, м3-мПа-1-с-1;
Та - временной период моделирования, с (сут);
V - элементарный объем, м3;
V - объем матрицы (породы), м3;
х, у - пространственные координаты, м; ф - пористость трещин, д.ед.;
// - вязкость жидкости, мПа-с; т - шаг по времени, с (сут);
V - вектор градиента;
О, - разностная сетка по времени;
0.х - разностная сетка по пространству;
Т -конечно-разностный оператор;
д - частная производная.
Верхний индекс:
0 - начальные или определенные условия; п - временной слой;
5 - номер итерации по нелинейности.
Нижний индекс:
1, ] - узлы пространственной разностной сет-
ки.
Список литературы
1. Quandalle P. Typical features of a multipurpose reservoir simulator / P.Quandalle, J.C.Sabathier // (SPE 16007) SPE Reservoir Engineering. November, 1989. P. 475-480.
2. Sabathier J.C. A new approach of fractured reservoirs / J.C.Sabathier, B.J.Bourbiaux., M.C.Cacas, S.Sarda // (SPE 39825) SPE International Petroleum Conference and Exhibition, Villahermosa, Mexico. 1998. P. 1-12.
3. Dershowitz W. Stochastic discrete fracture modelling of heterogeneous and fractured reservoirs / W.Dershowitz, N.Hurley, K.Been // Proceedings of the Third European Conference on the Mathematics of Oil Recovery. Delft, The Netherlands, 1992. P.119-135.
4. Dershowitz B. Integration of discrete feature network methods with conventional simulator approaches / B.Dershowitz,, P.LaPointe, T.Eiben, L.Wei // (SPE 62498) SPE Reservoir Evaluation & Engineering. April 2000. Vol.3. Num. 2. P.165-170.
5. Bourbiaux B.J. An integrated workflow to account for multi-scale fractures in reservoir simulation models: implementation and benefits / B.J.Bourbiaux, R.Basquet, M.C.Cacas, J.M.Daniel, S.Sarda // (SPE 78489) 10th Abu Dhabi International Petroleum Exhibition and Conference, 2002. P. 1-13.
6. Basque t R.. Fracture flow property identification: an optimized implementation of discrete fracture network models / R.Basquet, C.E.Cohen, B.Bourbiaux // (SPE 93748) 14th SPE Middle East Oil & Gas Show and Conference, Bahrain. 2005. P. 1-9.
7. Bogdanov I.I. Pressure drawdown well tests in fractured porous media / I.I.Bogdanov, V.V.Mourzenko, J.-F.Thover, P.M.Adler // Water resources research. Vol.38. 2002. P.X-1-X-19.
8. Викторин В.Д. Влияние особенностей карбонатных коллекторов на эффективность разработки нефтяных залежей /
В.Д.Викторин. М.: Недра, 1988. 150 с.
9. Shchipanov A.A. Simulation of fluid flow in fractured reservoirs subject to deformability / A.A.Shchipanov, A.Iu.Nazarov // AAPG International Conference & Exhibition. Proceedings. Paris, France, 2005. P. 1-6.
10.Голф-Рахт Т.Д. Основы нефтепромысловой геологии и разработки трещиноватых коллекторов: пер. с англ. / Т.Д.Голф-Рахт. М.: Недра, 1986. 608 с.
11.Азиз Х.Математическое моделирование пластовых систем: пер. с англ. / Х.Азиз,
Э.Сеттари. М.: Недра, 1982. 407 с.
12. Warren J.E. The behaviour of naturally fractured reservoirs / J.E.Warren & P.J.Root // SPE Journal. 1963. September. P.245-255.
A model of fluid flow in compressible discrete fracture network
A. A. Shchipanov, S. V. Rusakov
Perm State University, 614990, Perm, Bukireva st., 15
A new model of compressible discrete fracture network (CDFN) has been proposed to simulate fluid flow in deformable fractured rocks. Flow of homogeneous fluid in network of vertical intersecting fractures is considered. The fractures split impermeable rock into blocks that constitute matrix. Bond of fracture aperture and permeability is introduced from juxtaposition of Poisel's solution and Darcy's law. The matrix (rock) compressibility generates changing of fracture aperture due to pressure evolution during fluid flow. Proposed mathematical model is a basis for numerical simulation of flow in CDFN. In two-dimensional case a finite difference scheme and a research tool have been developed. The scheme is built based on finite volume method. The tool allows to simulate fluid flow in CDFN to a source (production or injection well). The model and the research tool can be used to characterize deformable fractured reservoirs of hydrocarbons, in particular to analyze transient pressure well test.