CC BY
111
2018, Т. 160, кн. 1 С.165-182
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
УДК 519.635.6
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЗАДАЧИ ДВУХФАЗНОЙ ФИЛЬТРАЦИИ В НЕОДНОРОДНЫХ ТРЕЩИНОВАТО-ПОРИСТЫХ СРЕДАХ С ИСПОЛЬЗОВАНИЕМ МОДЕЛИ ДВОЙНОЙ ПОРИСТОСТИ И МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ
В.И. Васильев, М.В. Васильева, А.В. Григорьев, Г.А. Прокопьев
Северо-Восточный федеральный университет имени М.К. Аммосова, г. Якутск, 677000, Россия
Аннотация
В работе проводится численное моделирование задач двухфазной фильтрации в трещиновато-пористых средах с использованием модели двойной пористости с сильно неоднородным коэффициентом проницаемости. Приводится система уравнений для случая двухфазной фильтрации без учета капиллярных и гравитационных сил, которая представляет собой связанную систему уравнений для давления и насыщенности в пористой среде имеющей систему трещин. Рассматриваются различные варианты задания функций перетока между пористой средой и трещинами. Численная реализация для аппроксимации скорости и давления строится на основе метода конечных элементов. Для дискретизации уравнения насыщенности посредством метода введения искусственной диффузии используется классический метод Галеркина с противопотоковой аппроксимацией. Приводятся результаты численных расчетов для модельной задачи с использованием различных функций перетока.
Ключевые слова: двухфазная фильтрация, неоднородные среды, трещиновато-пористые среды, модель двойной пористости, функции перетока, метод конечных элементов, численная стабилизация
Введение
В настоящее время запасов нефти в традиционных песчаных коллекторах практически не осталось, что приводит к необходимости разработки нетрадиционных коллекторов, которые характеризуются сложным неоднородным строением и наличием сети трещин (так называемые естественно-трещиноватые коллекторы). В подобных коллекторах матрица породы в большинстве случаев отличается низкой проницаемостью, поэтому течение происходит в основном в системе трещин. Отметим также, что гидродинамические характеристики реальных пластов существенно неоднородны.
При построении математической модели и проведении численного моделирования система трещин требует специального подхода, поскольку фильтрация в таких средах обладает специфическими особенностями. Трещины отличаются большой проницаемостью и оказывают существенное влияние на процессы переноса и течения в пористой среде. При этом следует разделять системы трещин, так как они могут сосуществовать на различных масштабах (микротрещины с размерами в несколько сантиметров, макротрещины, разломы), могут различаться природой
их возникновения (естественно-трещиноватые среды, разломы, трещины, возникшие за счет использования технологии гидроразрыва пласта). В случае естественно-трещиноватых пористых сред система трещин является в основном связанной, и для ее моделирования традиционно используют модели двойной пористости и мультиконтинуума [1—5]. Подобные модели строятся для идеализированной среды и обладают рядом ограничений. Взаимодействие континуумов в таких моделях задается посредством введения функций перетока между матрицей породы и трещинами. Определение этих функций является ключевой проблемой моделей мультиконтинуума, поскольку выбор обменных членов основывается на дополнительных допущениях.
В общем случае матрица пористой среды может быть неоднородной [6-8], что ведет к значительным вычислительным сложностям и, как было отмечено выше, обменный переток между средами только ухудшает общую картину. Кроме того, система трещин может содержать участки как с низкой, так и с высокой проницаемостью [9].
Основной идеей, на которой базируются модели двойной пористости, является раздельный расчет течения жидкостей в матрице и в трещинах с заданием обменного перетока между ними посредством использования функции переноса [10]. Данная модель была впервые предложена Баренблаттом [1]. Уоррен и Рут предложили модель [11], в которой все фильтрационные потоки идут в системе трещин, а поровая матрица участвует в этом процессе вторично через обмен. Далее Ка-земи [12] расширил модель Уоррена-Рута для двухфазной жидкости. Позднее в работах Дугласа и Арбогаста [3, 13] были рассмотрены модели однофазной и многофазной фильтрации для двойной пористости, построенные в результате использования процедуры усреднения.
В случае же несвязанных трещин, возникающих, например, при гидроразрыве пласта, применяемом для слабопроницаемых коллекторов, а также для больших относительно размеров пласта, несвязанных трещин или разломов следует строить согласованные сетки и использовать модели прямого учета течения в трещинах, такие как дискретная модель трещин (Discrete Fracture Model, DFM) или встроенная модель трещин (Embedded Fracture Model, EFM). Дискретная модель трещин позволяет учитывать особенность, связанную с трещинами, на уровне сетки и задавать на ней соответствующие свойства [14-16]. Данный подход является вычислительно трудоемким, поскольку подобный учет трещин при аппроксимации задачи приводит к большому количеству неизвестных. Следует также отметить, что некоторые трещины, размерности которых очень малы по сравнению с размерами пласта, невозможно точно захватить в узлы сетки. В отличие от метода DFM, метод EFM позволяет избежать условия согласованности сеток, и уменьшить размерность решаемой системы [17-20].
В настоящей работе рассматривается процесс двухфазной фильтрации в неоднородных трещиновато-пористых средах. Для численного моделирования используется модель двойной пористости с различными вариантами задания функций перетока. В работе используются классический и обобщенный виды перетоков, а также предлагается альтернативный (смесевыравнивающий) вариант. Смесевырав-нивающий альтернативный переток базируется на концепции баланса жидкостей в смеси. Если одна из жидкостей находится в переизбытке, то смесь стремится к балансу путем уменьшения доли данной жидкости и увеличения доли другой. В первой части работы представлена математическая модель, описывающая течения двухфазной жидкости в неоднородных трещиновато-пористых средах с использованием модели двойной пористости без учета гравитационных и капиллярных сил. Во второй части представлены различные варианты функций перетока: классический, обобщенный и альтернативный. Численная аппроксимация задачи
строится с помощью метода конечных элементов и представлена в третьей части статьи, где для уравнений переноса используется численная стабилизация посредством метода SUPG (streamline-upwind Petrov/Galerkin). Результаты численных расчетов для предложенных вариантов функции перетока представлены в четвертой части работы.
1. Математическая модель
Математическая модель течения жидкости в пористых средах базируется на законе сохранения массы для каждой из фаз
-—--+ V • (ри ) - Р1(11т = р^дз, а =
и обобщенном законе Дарси
иа = -ка gradра,
Рг
где значения индекса а = ], т относятся соответственно к трещинам и матрице пористой среды.
Здесь фа - пористости, иа - скорость течения г-й фазы, рг, р", Б", рг, кгг -плотность, давление, насыщенность, вязкость и относительная проницаемость г -й фазы, дз - дебиты добывающих/закачивающих скважин, к - проницаемость пористой среды, д"т - функции перетока между матрицей пористой среды и трещинами, ^^ Б" = 1. Отметим, что в статье мы пренебрегаем гравитационными
и капиллярными силами.
Рассмотрим далее случай двухфазной фильтрации нефти и воды (i = o,w) при фа = const и pi = const. При пренебрежении капиллярными силами давление для различных фаз является одинаковым
а а а и
Po = Pw = p , а = f,m-Пусть иа - суммарная скорость двухфазного потока
иа = иа + uW = -X(Sf)ka grad ра,
где Аа = \а + XW - общая подвижность, Xf - подвижность i -й фазы и ff - доля а -й фазы в двухфазном потоке
fa = А±_ ха = kri(Si )
Аа Pi
и иа = fWиа .
Таким образом, получаем следующую систему уравнений
d_sw dt
d_Sl dt
фа^Г + V^ (fw (S а)иа) - qWm = qW,
Фа^о + v^ (fo(sа)иа) - gfm = qf-
Поскольку ff + ff = 1, то складывая эти два уравнения, получаем
bsw
~3t
д sw
Фа — + V^ (fw (S а)иа) - qfm = qf, div иа - qff = qf
Рис. 1. Пористая среда и сеть трещин
где да = яОО + яОО, яО = Я0т + Я0т, в качестве расчетной насыщенности будем использовать насыщенность водной фазы Ба = БО.
Будем считать, что линейные размеры трещин достаточно малы по сравнению с характерными размерами пласта, в этом случае для их учета можно применить модель двойной пористости (рис. 1). Для моделей двойной пористости мы получим следующие системы:
для матрицы пористой среды
ЯСт
фт^)Г + У ' (Бт)ит) - Ягт = 0,
ит - дт = 0, ит + \(Бт)кт gradрт =0,
для сети трещин
* дБ1 * ,
ФТ + ^ ' (1ш (Б1 )и' ) + Ягт = Яг,
div и1 + Ят = Я,
и1 + \(Б( )к1 grad р1 = 0,
где
яО + Я°т = 0 Я°т = Ягт, яО = -Ятт = -Ят г =
Отметим, что в данной модели будем предполагать, что источниковые члены правой части задаются для сети трещин, что влияет в дальнейшем на поровую матрицу посредством обменного перетока.
2. Функции перетока
Рассмотрим функции перетока. В настоящей работе мы будем использовать следующие классические функции перетока, предложенные в работе Каземи [12]:
Ягт = (Бт)(рт - р1), Я от = я\0(Бт)(рт - р1).
Здесь а - коэффициент перетока между матрицей пористой среды и трещинами
а = + | + I) »л» а = 4(| + | + ^
где кш = (кш, 0, 0; 0, кш, 0; 0,0, кш) - абсолютная проницаемость пористой среды и 1х, 1у, - пространственные параметры трещин по направлениям х,у и я. В качестве функции перетока рассмотрим также обобщенный случай
Чшш = ап-ш (Рт - Р? ), Чаш = аЩо(рш - Р?), (1)
где щ характеризует направленный переток, щ = щг2 :
(ад^), И (р? -рш) > 0
Пг2 = { (2)
[Аг(Бш), И (р? - рш) < 0.
В качестве альтернативного перетока рассмотрим смесевыравнивающий переток. Такой переток позволяет учитывать процессы, связанные с обменом жидкости по составу происходящие в континуумах (подсредах), в результате которого происходит не только выравнивание давлений, но и выравнивание состава смесей,
Пг = Пг3 :
Пгз = sign (р? - рш)(Аг(Б?) - Аг(Бш)).
Второй альтернативный переток основан на среднеарифметическом значении проницаемости в трещинах и в матрице
Пг4 = (Аг(Б Ш)+ Аг(Б^ ))/2.
Отметим что для первого случая щ = щ1 = Аг (Бш). Поскольку
V • (и (Ба)иа) = паУ/ш (Ба) + и (Ба)У ■ иа, щ = щш + По, и? = ч - Чш, ^у иш = ч + Чш, Чш = Чшш + Чаш = ащ(рш - Р?), Чшш = ащш (Рш - Р?) = щЧш,
получим следующие уравнения для насыщенностей для матрицы пористой среды и для сети трещин:
ЯБш
Фш+ и^/ш(Бш) + и(Бш)чш - Чшш = 0, дБ?
Ф?+ и? V(Б ?) - /- (Б? )Чш + Чшш = Чш - /ш (Б? )ч.
Таким образом, для матрицы пористой среды получаем следующую систему дифференциальных уравнений неэволюционного типа:
дБ ш
Фш д— + и^/ш (Бш)+ а ( /ш (Бш) )(рш - р?) =0, х е П,
г, . I у J ш \ ^ у 1 \ ./ ш V ^ у
СЛ \ П
ё1у иш - ащ(рш - р? ) = 0, х е П, (3)
иш + кшА(Бш) gradрш = 0, х е П,
и для сети трещин
Ф1 дБ + и1 (Б1) - (Б1) - ^ (рт - р1) = Яг - и (Б1 )я, х е П,
div и1 + *п(рт - р1 ) = /, х е П, (4)
и1 + к1 Х(Б1) gradр1 =0, х е П,
где / = Яг - /г(Б1 )я.
Система уравнений (3)-(4) дополняется начальным условием для насыщенности
Бт(х, 0)= Б0, Б1 (х, 0)= БО, х е П, (5)
и граничными условиями для давления
рт = рО, р1 = р1, х е гl,
рт = Рт, р1 = р), х е Г2, (6)
ит = и1 =0, х е дП/Г} /г2,
и насыщенности
Бт = БО, Б1 = Б}, х е Г1, Бт = БО, Б} = Б), х е Г2. (7)
Отметим, что для простоты изложения в дальнейшем будем считать, что Я = Яо = = Яг = 0.
3. Конечно-элементная аппроксимация и вычислительный алгоритм
Для численного решения полученной начально-краевой задачи для системы уравнений неэволюционного типа (3)—(7) построим численную аппроксимацию задачи с помощью метода конечных элементов. Для дискретизации компонент скорости и давления по пространственным переменным применим смешанный метод конечных элементов с использованием элементов Тома - Равиара наименьшего порядка для скорости и кусочно-постоянных для давления. Данный метод позволит нам наиболее естественным образом аппроксимировать функции перетока, поскольку давление является в данном случае кусочно-постоянным. Смешанный метод конечных элементов дает нам также локально консервативную аппроксимацию, что является одним из необходимых условий при моделировании задач фильтрации.
Для аппроксимации рассматриваемой системы уравнений введем следующие пространства для скорости, давления и насыщенности:
ш = {т е ь2(п)а ^V е ь2(п)}, д = Ь2(П), V = Н2(П).
Вариационная постановка задачи записывается следующим образом для пористой среды
/ д Б т \ Л
™т{-д^, V0) + Ст(Бт, Vт) + яО(рт - р1, Vт) =0 У Vт е У,
8т(ит,тт) + Ь(тт,рт)=0 Утт е Ш, Ь(ит,ят) + я(рт -р1 ,чт) = 0 У ят е д,
и для сети трещин
т*(^У) + с^У) - ^ -* у)=о ^ е V,
в * (и*) + Ь(р*) = 0 V т* е Ш, Ь(и*, д*) - д(рт - р* ,д* )=0 V д* е ¿V,
где
уа) = ! фаБауа ¿х, са(Ба, уа) = ! ¡'а(Ба)иа ■ УБауа ¿х, п п
Ь(ра, иа ) = иара ¿х, ва(иа,та) = ! (ка\(Ба))-1иата ¿х,
д(рт - р* ,да) = ! ап(рт - р* )да ¿х, п
д* (рт - р*,Уа) = I а ¡г (Ба) + ^ (рт - р*)уа ¿х. п
где 51а е V, иа е Ш, ра е Я и а = т, /.
Аппроксимация уравнения для насыщенности с использованием стандартного метода Галеркина приводит к возникновению осцилляций в решении задачи. При аппроксимации уравнений конвекции-диффузии с использованием метода конечных разностей или метода конечных объемов в случае задач с преобладающим конвективным слагаемым, обычно применяются противопо-токовые аппроксимации конвективного слагаемого. Отметим также, что подобные схемы, в отличие от схем с центральными разностями, имеют первый порядок точности по пространственным переменным. Противопотоковые схемы можно построить посредством добавления искусственной диффузии для центрально-разностной схемы. Для конечно-элементной аппроксимации мы будем строить аналогичные противопотоковые схемы посредством введения дополнительного слагаемого (искусственной диффузии).
Пусть - триангуляция области П на элементы К с диаметрами Нк и к = = тах Нк. Для дискретизации по времени мы применим неявную аппроксимацию
по времени с шагом т и Ьп = пт. Запишем дискретную задачу следующим образом: для пористой среды
, от _ от )
тт[к т к ,ут) +Ст(Бт,ут) + ат(Бт,Ут) + дт (рт - р[,ут)=0 V Ут е я, вт(ит,тт) + Ь(рт,тт)=0 Vтт е (8)
Ь(ит,дт) + д(рт -р[,дт) = 0 Vдт е ¿¿н,
для сети трещин
т/( 0 - 0 +с* (Б[ у ) + а* (Б[ у) - д* (рт - р*У )=0 V V* е
в*(и{,т*)+ Ь(р{,т*) = 0 Vт* е Ь(и{,д*) - д(рт - р{,д*) = 0 V д* е
и
где п - номер временного слоя, Ба и Ба - насыщенность для текущего п-го слоя и предыдущего (п - 1)-го слоя соответственно. Здесь
аа(Ба, уа) = ! О^Ба ■ Vva ¿х,
где Оа = \иа\Н/2, \иа\ = -у/(иа, иа) - дополнительное слагаемое (искусственная диффузия), которое пропорционально шагу по сетке и скорости фильтрации [6, 21].
Систему уравнений (8)-(9) запишем в матричной форме (1 Мт + Ст(Бш,иш) + Бш(иш^ Бш + Я3ш(Бш, щш,По)(рш - Р?) = 1 МтБш,
,Пш ,По)(Р - Р'
Аш(Бт )ит + Врш = 0,
ВТиш + Q(пo,пш)(рш -Р?) = 0,
(10)
М? + С? (Б? ,и?)+ Б? (и? ^ Б? - Qsf (Б? ,Щш ,По)(рш - Р?) = 1 М? Б?,
А? (Б ? )и? + Вр? =0,
ВТи? - Q(пo,пш)(рш -Р?) = 0,
где щг = Пг(Б?,Бш,р?,рш) для г = о, т.
Полученная система уравнений (8)-(9) является сильно связанной, поскольку 1) перенос и течение в трещинах и пористой среде связаны за счет функций перетока, 2) уравнение для насыщенности зависит от скорости фильтрации и 3) уравнение для скорости и давления зависят от насыщенности [6, 10]. В случае линейных коэффициентов
Аш(Ба) = Ба, Ао(Ба) = 1 - Ба, А(Ба) = 1
связанную систему уравнений (10) удается расщепить по переменным. Тогда получим следующий вычислительных алгоритм.
• Вычислить скорости и давления для пористой среды и трещин
0 иш
( Аш В 0
ВТ Q 0
0 0 А?
0 ВТ
р"
и? \р?)
0,
где матрицы Аа, В и Q не зависят от насыщенности и не меняются по времени.
Вычислить насыщенности для каждого временного слоя Ьп < Ттах, п = = 1, 2,...
(1 Мш + Сш(Бш, иш) + Бш(иш)) Бш + QSш(Бш, Щш)(рш - Р?) = 1 МтБш, М? + С? (Б? ,и?)+ Б? (и? ^ Б? - QSf (Б? ,Щш)(рш - Р?) = 1 М? Б?.
Отметим, что в данном случае уравнение для насыщенности является нелинейным. Для его численного решения можно использовать те или иные методы линеаризации, например метод Ньютона, или в простейшем случае взять значения коэффициентов с предыдущего временного слоя (метод замороженных коэффициентов).
В случае нелинейных коэффициентов и Ха расщепить переменные скорости и давления с переменной насыщенности не представляется возможным. Поэтому мы использовали метод последовательного итерационного уточнения решения SEQ (sequantional), в котором последовательно находятся переменные (иа,ра) и Ба, вычисление которых производятся в рамках нелинейных связывающих итераций. Таким образом, вычислительный алгоритм можно записать следующим образом.
• Ьля каждого временного слоя Ьп < Ттах, п = 1, 2,...
— для каждой итерации к = 0,1,...
* вычислить скорости и давления (иа'к+1, ра'к+г), а = /, т:
[Ат(Бт°) В 0 0 \ /ит\
вТ ,пкк,) 0 ^(пк )
0 0 А1 (Б ?'к) В
\ 0 ^(пк ,пк) ВТ Q(r|k ,пк))
где пк = пк (Бт'к,Б,рт'к,рт'к) для % = о,т;
* вычислить насыщенности Ба = Ба'к+1
Ыт+Ст(Бт,ит'к+1) + Вт(ит'к+1)^ Бт+
+ QÍ(Бm,vkw+1,vko+1)(pm'k+1 - Р1к+1) = 1 МтБт, ^ М1 +С1 (Б1 ,и1'к+1) + (и1к+1 ^ Б1 -
- QSf (Б 1 ,пк+1,пк+1)(рт'к+1 - р10+1) = ~ М}Б1.
т
где п0+1 = п0+1 (Бт,Б1 ,рт,к+1,рт,к+1) для % = о, ю.
Данный алгоритм позволяет расщепить вычисление насыщенности для пористой среды и трещин.
\pf )
4. Результаты численных расчетов
Приведем результаты численного моделирования модели двойной пористости для фильтрации двухфазной жидкости. Дискретизация задачи приводится для обезразмеренного случая. Рассмотрим случай, когда qo = qw = 0, pw = = 1,
kwr (Sa) = Sa, kor (Sa) = 1 - Sa Xw(Sa) = Sa, Xo(Sa) = l - Sa, X(Sa) = l,
следовательно, fW(Sa) = Sa, fW = 1.
Численное моделирование проводилось при следующих значениях гидродинамических характеристик: пористость трещиноватой среды фf = 0.2, пористой среды фт =0.7, проницаемость поровой матрицы km = 1, проницаемость трещиноватой среды была получена из 77-го слоя теста SPE10 (см. рис. 2). Тест SPE10 представляет собой бенчмарк для апробации симуляторов добычи нефти. Модель имеет
Рис. 2. Абсолютная проницаемость трещин, в логарифмической шкале
Рис. 3. Распределение давления с использованием обобщенного перетока для Варианта 1 (<71 = 20). Слева - трещины; справа - поровая матрица
простую геометрию без особенностей. Причиной такой геометрии является обеспечение максимальной гибкости при выборе масштабируемых сеток. На подробной геологической структуре модель описывается с помощью обычной декартовой сетки. Размеры модели - 1200 х 2200 х 170 (в футах). Верхние 70 футов (35 слоев) представляют собой образования Тарберта, а нижние 100 футов (50 слоев) - образования Верхнего Несса. Размер ячейки подробной сетки составляет 20 х 10 х 2 фута. В тесте SPE10 всего 85 слоев, вместе они образуют трехмерный тест для апробации симуляторов разработки нефтяных месторождений. Проницаемость задается в виде упорядоченного текстового файла, из которого можно считать интересующие нас данные.
Пусть О - прямоугольная область 1.5x0.5, а Г и Г2 - левая и правая границы, соответственно. На левой границе задавались следующие граничные условия:
Sf =0.8, Sm = 0.2, pf = 1, pm = 1 ж e Гь
на правой границе давление задавалось равным нулю как для матрицы пористой среды, так и для трещин, на верхней и нижней границах области О было задано условие непроницаемости, то есть скорость фильтрации задавалась равной нулю. Расчет проводился для Tmax = 0.05 с использованием 200 временных слоев с т = = 2.5 • 10~4 (tm = тт, где т = Tmax/M, M = 200 и m - номер временного слоя).
Для проведения численных расчетов была построена неструктурированная расчетная сетка содержащая 7701 вершин и 15000 элементов. Для функции перетока зададим о\ = 20 (Вариант 1) и 02 = 100 (Вариант 2).
Приведем результаты численного моделирования задачи с использованием обобщенного перетока, представленного формулами (1), (2). На рис. 3-4 показаны распределения давления и насыщенности на различные моменты времени для Варианта 1, где о\ = 20. На рис. 5-6 показаны распределения давления и насыщенности для Варианта 2 с 02 = 100. Поскольку проницаемость для сети трещин принималась неоднородной, фронт насыщенности для трещин движется по каналам с высокой проницаемостью, и в этих местах также возникает переток в матрицу пористой среды. При этом если сравнить результаты для Варианта 1
Рис. 4. Насыщенность в различные моменты времени = 0.001, 0.01 и 0.05 (сверху вниз) с использованием обобщенного перетока для Варианта 1 (а± = 20). Слева - трещины; справа - поровая матрица
Рис. 5. Распределение давления с использованием обобщенного перетока для Варианта 2 (<2 = 100). Слева - трещины; справа - поровая матрица
Рис. 6. Насыщенность в различные моменты времени = 0.001, 0.01 и 0.05 (сверху вниз) с использованием обобщенного перетока для Варианта 2 (<2 = 100). Слева - трещины; справа - поровая матрица
Рис. 7. Насыщенность в различные моменты времени = 0.0001, 0.0015 и 0.02 (сверху вниз) с использованием обобщенного перетока для нелинейного случая. Слева - трещины; справа- - поровая матрица
(<71 = 20) и для Варианта 2 (<2 = 100), то в Варианте 2 наблюдается больший переток в матрицу за счет большего значения параметра 7 в функции перетока, распределения насыщенностей для трещин также отличаются.
Для демонстрации возможности обобщения результатов на случай нелинейных коэффициентов приведем результаты численного моделирования задачи с использованием обобщенного перетока, задаваемого формулами (1), (2), но в данном случае речь идет о нелинейных коэффициентах Х^ = в2/^. Для реализации нелинейности был использован последовательный итерационный алгоритм Ньютона. Количество итераций метода Ньютона для алгоритма варьировалось от 2 до 5 для давления - скорости и от 4 до 10 для насыщенности. На рис. 7 представлены распределения насыщенности на различные моменты времени. Остальные коэффициенты представлены следующими значениями: вязкость воды = 0.0003, вязкость нефти 1ла = 0.003, mf = 0.2, тт = 0.7, < = 100. Для численного моделирования нелинейного случая был выбран слой 77 из теста 8РЕ10. Вычислительные параметры сетки и шага по времени: т = 0.0001, Ттах = 0.02, расчетная сетка состоит из 5243 вершин и 10196 элементов. Рис. 8-9 иллюстрируют интегральную характеристику по правой границе области для Варианта 1 и Варианта 2, вычисленной по формуле
= ¿У < • п ¿в, п1 = /апа при /а = а = т, /.
¿0 Р
г 2
Здесь суммирование производится по слоям по времени до момента Ьт, п - нормаль к границе Г2 . Приведенные результаты характеризуют обводненность добывающей скважины, использование большего значения < в Варианте 2 приводит к большим значениям Qa как для матрицы пористой среды, так и для сети трещин (см. рис. 8-9). Отметим, что интегральные характеристики для классического, обобщенного и усредненного перетоков практически совпадают, а альтернативный переток приводит к меньшей обводненности в трещинах, но к большей для матрицы пористой среды.
Рис. 8. Интегральная характеристика для насыщенности на правой границе Г2 для различных вариантов функций перетока для Варианта 1 (<1 = 20). Слева - трещины; справа - поровая матрица
Рис. 9. Интегральная характеристика для насыщенности на правой границе Г2 для различных вариантов функций перетока для Варианта 2 (<2 = 100). Слева - трещины; справа - поровая матрица
Исходя из полученных результатов расчетов, можно сделать следующие выводы:
- выполнена численная реализация математической модели двойной пористости для двухфазной жидкости в трещиновато-пористой гетерогенной среде;
- трещины служат каналами высокой проводимости, из которых потом флюид перетекает в матрицу;
- рассмотрены различные варианты перетоков, для всех перетоков модель продемонстрировала свою работоспособность.
Следует отметить, что в настоящей работе был рассмотрен некоторый ограниченный класс функций перетока, в будущем будут рассмотрены и другие варианты функций перетока, такие как градиентный, модель MINC (Multiple INteracting Continua).
Заключение
В работе рассмотрена задача двухфазной фильтрации жидкости для неоднородных трещиновато-пористых сред. Математическая модель описана системой связанных уравнений для насыщенностей, скоростей и давлений в матрице пористой среды и в сети трещин на основе модели двойной пористости. Для численного решения была построена конечно-элементная аппроксимация по пространству
с использованием элементов Тома - Равиара для переменных скоростей и давлений, и классического метода Галеркина со стабилизацией для насыщенностей. В дальнейшем мы планируем рассмотреть более общие модели с учетом капиллярных и гравитационных сил и построить численное решение задачи для трехмерных квазиреальных геометрий.
Благодарности. Работа выполнена при финансовой поддержке РФФИ (проект № 17-01-00732) (В.И. Васильев: теория, постановка задачи), Мегагранта правительства РФ N 14.Y26.31.0013 (М.В. Васильева: численный алгоритм; Г.А. Прокопьев: реализация нелинейного случая), а также РНФ (проект № 17-71-10106) (А.В. Григорьев: реализация линейного случая, эксперименты).
Литература
1. Баренблатт Г.И., Желтое Ю.П., Кочина И.Н. Об основных представлениях теории фильтрации однородных жидкостей в трещиноватых породах // Прикл. матем. и механика. - 1960. - Т. 24, Вып. 5. - С. 852-864.
2. Bourgeat A. Homogenized behavior of two-phase flows in naturally fractured reservoirs with uniform fractures distribution // Comp. Methods Appl. Mech. Eng. - 1984. - V. 47, No 1. - P. 205-216.
3. Arbogast T., Douglas Jr. J., Hornung U. Derivation of the double porosity model of single phase flow via homogenization theory // SIAM J. Math. Anal. - 1990. - V. 21, No 4. -P. 823-836.
4. Gwo J.P., Jardine P.M., Wilson G.V., Yeh G.T. A multiple-pore-region concept to modeling mass transfer in subsurface media // J. Hydrol. - 1995. - V. 164, No 1. -P. 217-237.
5. Kalinkin A.A., Laevsky Y.M. Mathematical model of water-oil displacement in fractured porous medium // Sib. Electron. Math. Rep. - 2015. - V. 12. - P. 743-751.
6. Васильева М.В., Васильев В.И., Тимофеева Т.С. Численное решение методом конечных элементов задач диффузионного и конвективного переноса в сильно гетерогенных пористых средах // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2016. -Т. 158, кн. 2. - C. 243-261.
7. Васильева М.В., Васильев В.И., Красников А.А., Никифоров Д.Я. Численное моделирование течения однофазной жидкости в трещиноватых пористых средах // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2017. - Т. 159, кн. 1. - C. 100-115.
8. Ginting V., Pereira F., Presho M., Wo S. Application of the two-stage Markov chain Monte Carlo method for characterization of fractured reservoirs using a surrogate flow model // Comput. Geosci. - 2011. - V. 15, No 4. - P. 691-707. - doi: 10.1007/s10596-011-9236-4.
9. Salimi H., Bruining J. Improved prediction of oil recovery from waterflooded fractured reservoirs using homogenization // SPE Reservoir Eval. Eng. - 2010. - V. 13, No 1. -P. 44-55. - doi: 10.2118/115361-PA.
10. Вабищевич П.Н., Григорьев А.В. Численное моделирование фильтрации флюида в анизотропной трещиновато-пористой среде // Сиб. журн. вычисл. матем. - 2016. -Т. 19, № 1. - С. 61-74.
11. Warren J.E., Root P.J. The behavior of naturally fractured reservoirs // Soc. Pet. Eng. J. - 1963. - V. 3, No 3. - P. 245-255. - doi: 10.2118/426-PA.
12. Kazemi H., Merrill L.S. Jr., Porterfield K.L., Zeman P.R. Numerical simulation of water-oil flow in naturally fractured reservoirs // Soc. Pet. Eng. J. - 1976. - V. 16, No 6. -P. 317-326. - doi: 10.2118/5719-PA.
13. Douglas J. Jr., Arbogast T., Paes-Leme P.J., Hensley J.L., Nunes N.P. Immiscible displacement in vertically fractured reservoirs // Transp. Porous Media. - 1993. - V. 12, No 1. - P. 73-106. - doi: 10.1007/BF00616363.
14. Karimi-Fard M., Durlofsky L.J., Aziz K. An efficient discrete fracture model applicable for general purpose reservoir simulators // SPE Reservoir Simul. Symp. - Soc. Petr. Eng., 2003. - Art. 79699, P. 1-11. - doi 10.2118/79699-MS.
15. Efendiev Y., Lee S., Li G., Yao J., Zhang N. Hierarchical multiscale modeling for flows in fractured media using generalized multiscale finite element method // Int. J. Geomath. -2015. - V. 6, No 2. - P. 141-162. - doi: 10.1007/s13137-015-0075-7.
16. Akkutlu I.Y., Efendiev Y., Vasilyeva M.V. Multiscale model reduction for shale gas transport in fractured media // Comput. Geosci. - 2016. - V. 20, No 5. - P. 953-973. -doi: 10.1007/s10596-016-9571-6.
17. Karimi-Fard M., Gong B., Durlofsky L.J. Generation of coarse-scale continuum flow models from detailed fracture characterizations // Water Resour. Res. - 2006. - V. 42, No 10. - Art. W10423, P. 1-13. - doi: 10.1029/2006WR005015.
18. Gong B., Karimi-Fard M., Durlofsky L.J. Upscaling discrete fracture characterizations to dual-porosity, dual-permeability models for efficient simulation of flow with strong gravitational effects // SPE J. - 2008. - V. 13, No 1. - P. 58-67. - doi: 10.2118/102491-PA.
19. Lee S.H., Jensen C.L., Lough M.F. Efficient finite-difference model for flow in a reservoir with multiple length-scale fractures // SPE J. - 2000. - V. 5, No 3. - P. 1-11. - doi: 10.2118/65095-PA.
20. Li L., Lee S.H. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media // SPE Reservoir Eval. Eng. -2008. - V. 11, No 4. - P. 750-758. - doi: 10.2118/103901-PA.
21. Donea J., Huerta A. Finite Element Methods for Flow Problems. - John Wiley & Sons, 2003. - xii, 350 p. - doi: 10.1002/0470013826.
Поступила в редакцию 18.12.17
Васильев Василий Иванович, доктор физико-математических наук, профессор, заведующий кафедры «Вычислительные технологии»
Северо-Восточный федеральный университет имени М.К. Аммосова
ул. Белинского, д. 58, г. Якутск, 677000, Россия E-mail: vasvasil@mail.ru
Васильева Мария Васильевна, кандидат физико-математических наук, доцент-исследователь кафедры «Вычислительные технологии»
Северо-Восточный федеральный университет имени М.К. Аммосова
ул. Белинского, д. 58, г. Якутск, 677000, Россия E-mail: vasilyevadotmdotv@gmail.com
Григорьев Александр Виссарионович, кандидат физико-математических наук, доцент-исследователь кафедры «Вычислительные технологии»
Северо-Восточный федеральный университет имени М.К. Аммосова
ул. Белинского, д. 58, г. Якутск, 677000, Россия E-mail: re5itsme@gmail.com
Прокопьев Григорий Анатольевич, аспирант кафедры «Вычислительные технологии»
Северо-Восточный федеральный университет имени М.К. Аммосова
ул. Белинского, д. 58, г. Якутск, 677000, Россия E-mail: dju92@mail.ru
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online) UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)
2018, vol. 160, no. 1, pp. 165-182
Mathematical Modeling of the Two-Phase Fluid Flow
in Inhomogeneous Fractured Porous Media using the Double Porosity Model and Finite Element Method
V.I. Vasilyev * , M.V. Vasilyeva ** , A.V. Grigorev *** , G.A. Prokopiev ****
M.K. Ammosov North-Eastern Federal University, Yakutsk, 677000 Russia E-mail: *vasvasil@mail.ru, ** vasilyevadotmdotv@gmail.com,
***re5itsme@gmail.com, ****reilroot@gmail.com
Received December 18, 2017 Abstract
Numerical simulation of the two-phase fluid flow in a fractured porous media using the double porosity model with a highly inhomogeneous permeability coefficient has been studied. A system of equations has been presented for the case of two-phase filtration without capillary and gravitational effects, which is a connected system of equations for pressure and saturation in a porous medium that contains a system of cracks. Different variants of specifying the flow functions between the porous medium and cracks have been considered. The numerical implementation for velocity and pressure approximation is based on the finite element method. To discretize the saturation equation, the classical Galerkin method with counter-flow approximation has been used. The results of numerical calculations for the model problem using various interflow functions have been presented.
Keywords: two-phase filtration, inhomogeneous media, fractured-porous media, double porosity model, interflow functions, finite element method, numerical stabilization
Acknowledgments. The study was supported by the Russian Foundation for Basic Research (project no. 17-01-00732) (V.I. Vasilyev: theory, problem statement), Megagrant of the Government of the Russian Federation (grant no. 14.Y26.31.0013) (M.V. Vasilyeva: numerical algorithm; G.A. Prokopiev: nonlinear case implementation), and by the Russian Science Foundation (project no. 17-71-10106) (A.V. Grigorev: linear case implementation, experiments).
Figure Captions
Fig. 1. Porous medium and fracture network.
Fig. 2. Absolute permeability of fractures, kf on the logarithmic scale.
Fig. 3. Pressure distribution using the generalized interflow for Variant 1 (a1 = 20). On the left - fractures; on the right - porous matrix.
Fig. 4. Saturation at different points of time tm = 0.001, 0.01, and 0.05 (from top to bottom) using the generalized interflow for Variant 1 (<j1 = 20). On the left - fractures; on the right - porous matrix.
Fig. 5. Pressure distribution using the generalized interflow for Variant 2 (a2 = 100). On the left - fractures; on the right - porous matrix.
Fig. 6. Saturation at different points of time tm = 0.001, 0.01, and 0.05 (from top to bottom) using the generalized interflow for Variant 2 (a2 = 100). On the left - fractures; on the right - porous matrix.
Fig. 7. Saturation at different points of time tm = 0.0001, 0.0015, and 0.02 (from top to bottom) using the generalized interflow for the nonlinear case. On the left - fractures; on the right - porous matrix.
Fig. 8. The integral characteristic for saturation at the right-hand boundary r2 for different variants of the interflow functions for Variant 1 (a1 = 20). On the left - fractures; on the right - porous matrix.
Fig. 9. The integral characteristic for saturation at the right-hand boundary r2 for different variants of the interflow functions for Variant 2 (a2 = 100). On the left - fractures; on the right - porous matrix.
References
1. Barenblatt G.I., Zheltov Iu.P., Kochina I.N. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]. J. Appl. Math. Mech., 1960, vol. 24, no. 5, pp. 1286-1303. doi: 10.1016/0021-8928(60)90107-6.
2. Bourgeat A. Homogenized behavior of two-phase flows in naturally fractured reservoirs with uniform fractures distribution. Comput. Methods Appl. Mech. Eng., 1984, vol. 47, no. 1, pp. 205-216.
3. Arbogast T., Douglas Jr. J., Hornung U. Derivation of the double porosity model of single phase flow via homogenization theory. SIAM J. Math. Anal., 1990, vol. 21, no. 4, pp. 823-836.
4. Gwo J.P., Jardine P.M., Wilson G.V., Yeh G.T. A multiple-pore-region concept to modeling mass transfer in subsurface media. J. Hydrol., 1995, vol. 164, no. 1, pp. 217-237.
5. Kalinkin A.A., Laevsky Y.M. Mathematical model of water-oil displacement in fractured porous medium. Sib. Electron. Math. Rep., 2015, vol. 12, pp. 743-751.
6. Vasilyeva M.V., Vasilyev V.I., Timofeeva T.S. Numerical solution of the convective and diffusive transport problems in a heterogeneous porous medium using finite element method. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2016, vol. 158, no. 2, pp. 243-261. (In Russian)
7. Vasilyeva M.V., Vasilyev V.I., Krasnikov A.A., Nikiforov D.Ya. Numerical simulation of single-phase fluid flow in fractured porous media. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2017, vol. 159, no. 1, pp. 100-115. (In Russian)
8. Ginting V., Pereira F., Presho M., Wo S. Application of the two-stage Markov chain Monte Carlo method for characterization of fractured reservoirs using a surrogate flow model. Comput. Geosci., 2011, vol. 15, no. 4, pp. 691-707. doi: 10.1007/s10596-011-9236-4.
9. Salimi H., Bruining J. Improved prediction of oil recovery from waterflooded fractured reservoirs using homogenization. SPE Reservoir Eval. Eng., 2010, vol. 13, no. 1, pp. 44-55. doi: 10.2118/115361-PA.
10. Vabishchevich P.N., Grigoriev A.V. Numerical modeling of fluid flow in anisotropic fractured porous media. Numer. Analys. Appl., 2016, vol. 9, no. 1, pp. 45-56. doi: 10.1134/S1995423916010055.
11. Warren J.E., Root P.J. The behavior of naturally fractured reservoirs. Soc. Pet. Eng. J., 1963, vol. 3, no. 3, pp. 245-255. doi: 10.2118/426-PA.
12. Kazemi H., Merrill L.S. Jr., Porterfield K.L., Zeman P.R. Numerical simulation of water-oil flow in naturally fractured reservoirs. Soc. Pet. Eng. J., 1976, vol. 16, no. 6, pp. 317326. doi: 10.2118/5719-PA.
13. Douglas J. Jr., Arbogast T., Paes-Leme P.J., Hensley J.L., Nunes N.P. Immiscible displacement in vertically fractured reservoirs. Transp. Porous Media, 1993, vol. 12, no. 1, pp. 73-106. doi: 10.1007/BF00616363.
14. Karimi-Fard M., Durlofsky L.J., Aziz K. An efficient discrete fracture model applicable for general purpose reservoir simulators. SPE Reservoir Simul. Symp. Soc. Petr. Eng., 2003, art. 79699, pp. 1-11. doi: 10.2118/79699-MS.
15. Efendiev Y., Lee S., Li G., Yao J., Zhang N. Hierarchical multiscale modeling for flows in fractured media using generalized multiscale finite element method. Int. J. Geomath., 2015, vol. 6, no. 2, pp. 141-162. doi: 10.1007/s13137-015-0075-7.
16. Akkutlu I.Y., Efendiev Y., Vasilyeva M.V. Multiscale model reduction for shale gas transport in fractured media. Comput. Geosci., 2016, vol. 20, no. 5, pp. 953-973. doi: 10.1007/s10596-016-9571-6.
17. Karimi-Fard M., Gong B., Durlofsky L.J. Generation of coarse-scale continuum flow models from detailed fracture characterizations. Water Resour. Res., 2006, vol. 42, no. 10, art. W10423, pp. 1-13. doi: 10.1029/2006WR005015.
18. Gong B., Karimi-Fard M., Durlofsky L.J. Upscaling discrete fracture characterizations to dual-porosity, dual-permeability models for efficient simulation of flow with strong gravitational effects. SPE J., 2008, vol. 13, no. 1, pp. 58-67. doi: 10.2118/102491-PA.
19. Lee S.H., Jensen C.L., Lough M.F. Efficient finite-difference model for flow in a reservoir with multiple length-scale fractures. SPE J., 2000, vol. 5, no. 3, pp. 1-11. doi: 10.2118/65095-PA.
20. Li L., Lee S.H. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Eval. Eng., 2008, vol. 11, no. 4, pp. 750-758. doi: 10.2118/103901-PA.
21. Donea J., Huerta A. Finite Element Methods for Flow Problems. John Wiley & Sons, 2003, xii, 350 p. doi: 10.1002/0470013826.
Для цитирования: Васильев В.И., Васильева М.В., Григорьев А.В., Проко-I пьев Г.А. Математическое моделирование задачи двухфазной фильтрации в неоднородных трещиновато-пористых средах с использованием модели двойной пористости \ и метода конечных элементов // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. -2018. - Т. 160, кн. 1. - С. 165-182.
For citation: 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. Uchenye Zapiski \ Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2018, vol. 160, no. 1, pp. 165-182. (In Russian)
