Сетевое и континуальное моделирование капиллярной неравновесности двухфазной фильтрации

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

Классическим способом моделирования двухфазных течений в пористых средах является закон Дарси:

q = — - VP, (1)

Ц где q — скорость потока, VP — градиент давления, ц — коэффициент динамической вязкости, К = kok(S), ко - абсолютная проницаемость, k(S) — относительная проницаемость, S — насыщенность смачивающей жидкостью.

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

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

^k =ФИ)-

где к — относительная фазовая проницаемость, S — насыщенность смачивающей жидкостью. t — время.

Однако процесс перераспределения флюида в поровом пространстве может происходить даже при постоянной насыщенности S = const. Этот факт можно учесть, включив внутренний параметр релаксации £ в аргументы к. Относительная фазовая проницаемость принимает вид k = k(S,£),                                             (3)

где k — относительная фазовая проницаемость, S — насыщенность смачивающей жидкостью, £ — внутренний параметр.

Требуется кинетическое уравнение, такое что £ релаксирует к равновесному значению для данной постоянной насыщенности S:

Таблица!

Различные подходы к моделированию двухфазной фильтрации в пористых средах

Континуальные модели (макромасштаб)	Неконтинуальные модели (модели порового масштаба) учитывают неравновесные эффекты
Классические континуальные модели: закон Дарси, модель Бакли - Леверетта - не учитывают неравновесные эффекты	Сетевые модели - простые неконтинуальные модели с высокой скоростью расчетов
Продвинутые континуальные модели: модель Хассанизаде и Грея, модель Баренблатта неравновесной эффективной насыщенности, Континуальная релаксационная модель [7] - учитывают неравновесные эффекты	Прямое численное моделирование (DNS): метод решеточных уравнений Больцмана (LBM), методы конечных объемов/элементов (FVM/FEM) -решают уравнения Навье - Стокса в поровом масштабе, вычислительно неэффективно

Ранее авторами была разработана сетевая модель [13], которая учитывает капиллярную неравновесность и использует новый метод распределения двух различных фаз в узлах. С использованием этой сетевой модели была смоделирована пропитка изолированного пористого блока, состоящего из двух областей с контрастными капиллярными свойствами. Внешняя область, содержащая более широкие капилляры, изначально была насыщена смачивающей жидкостью, а внутренняя (с более тонкими капиллярами) — несмачивающей. Результаты показали, что из-за капиллярного давления смачивающая жидкость вытесняет несмачивающую и проникает во внутреннюю область. Обнаружено, что насыщенность внутренней области смачивающей жидкостью сначала быстро возрастает, затем скорость роста уменьшается и насыщенность стремится к равновесному значению. Это согласуется с поведением параметра ф который релаксирует к равновесному значению. Была получена зависимость капиллярного давления от различных конечных насыщенностей смачивающей жидкости во внутренней области, которая качественно совпала с классическими литературными данными. Таким образом, была подтверждена валидность модели для моделирования пропитки.

Затем было смоделировано двухфазное вытеснение, при котором среда изначально была насыщена несмачивающей жидкостью, а для её вытеснения закачивалась смачивающая жидкость. Среда была неоднородной с периодическим изменением пористости и проницаемости. Были получены графики зависимости насыщенности от координаты S (ж) в различные моменты времени для различных коэффициентов поверхностного натяжения а. Для течений с очень низким а флюиды в основном текли по пути наименьшего сопротивления (более толстые трубки, которые являются областями с высокой проницаемостью). Однако при высоком а течение происходило через области с низкой проницаемостью из-за эффекта пропитки, вызванного капиллярными силами.

В данной статье изучается двухфазное вытеснение в среде с постоянной пористостью и периодически переменной проницаемостью. Полученные с помощью сетевой модели результаты S(x), сравниваются с численным решением континуальной модели [7] для различных параметров, чтобы проверить валидность сетевой модели и понять определённые параметры в континуальной модели [7].

В разделе 2 описана сетевая модель и начальные условия двухфазного вытеснения.

Раздел 3 содержит описание континуальной модели [7], системы её дифференциальных уравнений и метод их решения. В разделе 4 представлены и сравниваются результаты сетевой и континуальной моделей. Выводы приведены в разделе 5.

2.    Двухфазное вытеснение с использованием сетевой модели

Сетевая модель представляет собой упрощение реальной пористой среды, содержащей узлы, соединённые друг с другом трубками. Базовая схема сетевой модели [13] приведена па рис. 1, здесь узел ио соединён с четырьмя другими узлами (щ, ^2, U3, U4) посредством трубок. Трубки имеют цилиндрическую форму и могут быть разного радиуса, а также могут содержать мениски. Жидкости предполагаются несжимаемыми, следовательно, в каждом узле выполняется закон сохранения объёма:

^Q i = 0, i = (1,2, 3,4),                                  (5)

i где Qi — объёмный расход в узел через трубку i.

Рис. 1. Схема сетевой модели: (1) - смачивающая жидкость (тёмный оттенок), (2) - несмачивающая жидкость (светлый оттенок), (3) - мениск, Qk - объёмный расход в трубке i, Ui - узел i

Расход в каждой трубке зависит от разности давлений между соединяемыми узлами, количества и ориентации менисков:

^Q jk ^a jk ^^P jk + ^b jk ,                                  (б)

где Q jk — объёмный расход из узла j в узел k, A jk и B jk — параметры трубки, соединяющей узлы j и k, a AP jk = P j — P k — разность давлений между узлами, a jk учитывает геометрию трубки, в то время как b jk — разность давлений, обусловленную наличием менисков. Более подробно выражения для коэффициентов уравнения (6) представлены в [13].

Пусть в сети имеется п узлов и P j обозначает давление в узле j, тогда j — целое число такое, что: 1 <  j ^ п. Сетевая модель обновляет положения флюидов после каждого временного шага. Для каждого временного шага выполняется следующее:

• 1)    Для определения давлений в каждом узле инициализируется структура данных следующего типа:

А ? = В ,                               (7)

здесь А — матрица размером п х п, р = (Pi, P2, ...,P_n ) — вектор давлений, а В — матрица размером 1 х п. Изначально матрицы заполняются нулями: А = 0, В = 0.

• 2)    Пусть M i,j обозначает элемент в i-й строке и j-м столбце матрицы М. Для каждого узла, j устанавливается j-я строка, матриц А и В.

• а)    Если j является открытым узлом, то давление P j в этом узле известно из гранич-HBix условий. Например, пуств P b _oun d ^- известное дявление, тогда P j = P b _oun d, что задаётся следующим образом:

A j,j ^— 1, B j, 1 = P bound .

• б)    Если j является закрытым узлом, то к этому узлу применяется закон сохранения объёма — Ур. (5). Например, пуств узел j соединён с узлами ink, тогда:

A jd —   ^a j-

A j,j ^{— a} j- + ^a jk ,

^A j,k      ^a jk ,

^B j, i ^{— —b} j- ^- b jk .

• 3)    Система линейных алгебраических уравнений (7) решается методом Гаусса, что поз

воляет определить давления в каждом узле.

• 4)    По известным давлениям в каждом узле определяется расход в каждой трубке с использованием Ур. (6).

• 5)    Скорость движения менисков и флюидов в трубке определяется по формуле

^Vi = ^R ,                             ^и0)

где v_i — средняя скорость псрсэющепия флюида в трубке i. Q_i — объёмный расход согласно ур. (6). R i — радиус трубки.

• 6)    Выбирается подходящий шаг интегрирования по времени:

At — c_t ( VV- ^    ,                                   (11)

где At — шаг по времени, c t — 0.1 используется для вынш/.тений в данной статье. I- — длина трубки i. V- — скорость течеппя флюида в трубке i.

• 7)    Интегрирование выполняется путём смещения менисков и флюидов в каждой трубке на величину Ax - — Atv-. Используется соответствующий масштаб, и уравнения применяются в безразмерной форме в соответствии с [13].

• 8)    Для случаев, когда две фазы одновременно поступают в узел, данная сетевая модель использует новый алгоритм распределения флюидов, согласно которому смачивающая жидкость распределяется первой по трубкам в порядке возрастания их радиусов. Пример нового метода показан на рис. 2. В течение временного шага были вычислены объёмные расходы Q - для каждой трубки i — (1, 2, 3, 4) и выбран соответствующий шаг по времени, показанный на рис. 2а. В узел из трубки 4 и трубки 3 поступает V_w — 0.7 и V_nw — 0.3 безразмерных единиц объёма (далее - ед. объёма) смачивающей и несмачивающей жидкости соответственно.

Q2At — 0.2 и QiAt — 0.8 ед. объёма некоторых флюидов вытекают из узла через трубку 2 и трубку 1 соответственно. Поскольку трубка 2 тоньше, в неё сначала распределяется смачивающая жидкость. Общее количество смачивающей жидкости, поступающей в узел, составляет V_w — 0.7 > Q2At, следовательно, трубка 2 получает только смачивающую жидкость, как показано на рис. 26. Оставшиеся ед. объёма смачивающей жидкости ( V_w — Q2At — 0.5) добавляются в трубку 1. В конце несмачивающая жидкость V_nw — 0.5 распределяется в трубку 1, так как это единственная оставшаяся трубка, которая может принять жидкость.

Рис. 2. Новый метод распределения фаз при одновременном поступлении двух различных фаз в узел: (1) - смачивающая жидкость, (2) - иесмачивающая жидкость, Qi - объёмный расход из трубки i, At - шаг интегрирования по времени, w - смачивающая жидкость, nw - иесмачивающая жидкость, V - объём жидкости, добавленный или удалённый, (а) - сначала определяется поток жидкостей в узел и из узла, (б) - затем смачивающая жидкость распределяется в более топкие трубки, и, наконец, распределяется иесмачивающая жидкость. Q, V и t являются безразмерными величинами

• 9)    Если в трубке находится более 2 менисков, то фазы объединяются таким образом, что центр масс каждой фазы остаётся неизменным. Это предотвращает превышение максимального количества менисков в трубке сверх 2.

• 10)    В новом распределении, полученном после перемещения флюидов, проводятся различные измерения, например: капиллярное давление, насыщенность смачивающей жидкостью в области.

С использованием данной сетевой модели моделируется двухфазное вытеснение. Смачивающая жидкость закачивается в двумерный блок пористой среды через открытую границу слева, в то время как верхняя и нижняя границы являются закрытыми, как показано на рис. За. Пористый блок изначально насыщен несмачивающей жидкостью, а флюиды удаляются через правую открытую границу. В каждый момент времени t измеряется S (ж) - насыщенность смачивающей жидкостью в зависимости от координаты ж. Координата ж направлена вдоль приложенного градиента давления и лежит на горизонтальной оси. Используются безразмерные координаты, такие что: 0 С ж С 1 и 0 С У С 1- Дискретное распределение S (ж) получается путём вычисления насыщенности смачивающей жидкостью в тонкой вертикальной полосе шириной Аж. Для упрощения вычислений эта вертикальная полоса содержит один столбец трубок, центры которых имеют одинаковую координату ж. Использовалась сеть размером 40 х 40 трубок, следовательно, Аж = 1/40.

Смачивающая жидкость закачивается под высоким давлением P i _n(t), которое вычисляется на каждом временном шаге для поддержания постоянного объёмного расхода закачки Q', процедура расчета приведена в [12]. На правой открытой границе поддерживается постоянное давление P_ou t = 0. Изменение проницаемости достигается за счёт периодического изменения радиуса трубок:

^(ж,у) = А

[ 1 + В cos

(тж) с“(тУ)]’

где А = 0.5, В = 0.5, А = 0.5, хну- центры трубок в плоскости ху. График этого распределения показан на рис. 36. Цвет каждой ячейки обозначает радиус трубки в данном месте. Области более тёмного цвета соответствуют областям с ввюокой проницаемостью, а более светлые - с низкой.

Рис. 3. Схема двухфазного вытеснения и распределения радиусов: (а) - смачивающая жидкость закачивается через левую открытую границу, Q - объёмный расход, P_in(t) - давление на входе, которое вычисляется для каждого шага интегрирования в момент времени t, P_out - давление на выходе, Ах - тонкая полоса, в которой измеряется насыщенность, R) ~ R ~ радиус трубки, R_min^_ минимальный радиус в системе

В каждый момент времени вычисляется распределение насыщенности S (х) по координате х (вдоль направления приложенного градиента давления). Это делается путём нахождения насыщенности смачивающей жидкостью в вертикальной полосе шириной Ах, также показанной на рис. За. В данной статье пористость постоянна по координатам за счёт использования трубок постоянного объёма:

сг

^li     R2 ■

^с г — ^V tube^/’ ^,

где l_i - длина трубки i, Ri - радиус трубки i, с_г и V_tu b _e - константы.

3.    Неравновесная континуальная модель

В задаче Бакли - Леверетта одномерная пористая среда изначально насыщена несмачивающей жидкостью, а смачивающая жидкость закачивается с левого конца. Требуется найти S(х■t) - насыщенность смачивающей жидкостью по координате х в каждый момент времени t. Уравнение переноса насыщенности в безразмерной форме имеет вид:

dS db(S) dS =₀

dt dS дх ■ где S = S(х■t') - насыщенность смачивающей жидкостью, t - безразмерное время, такое что при t = 1 в (•истому будет закачан один поровый объём смачивающей жидкости, х -безразмерная координата, такая что левый конец одномерной пористой среды находится при х = 0, а правый - при х = 1, b(S) - функция Бакли, такая что

ЫУ

b(S) = ЦМ^МУТ ■ где ki - относительная фазовая проницаемость смачивающей жидкости, к2 - относительная фазовая проницаемость несмачивающей жидкости, Ц1 - вязкость смачивающей жидкости, Ц2 - вязкость несмачивающей жидкости. Для расчетов в данной статье использованы зависимости:

ki = S 2; к2 = (1 -S)². (16)

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

S(xo , t) = 1;    S(х, to) = 0,

где хо = 0 11 to = 0.

Модель Бакли - Леверетта основана на законе Дарси и предполагает, что течение происходит очень медленно по сравнению с перераспределением флюида в поровом пространстве, которое заключается в вытеснении смачивающей жидкостью несмачивающей из более мелких пор под действием капиллярного давления. Охарактеризуем это время с помощью т, которое является безразмерной переменной с тем же масштабом, что и наше безразмерное время t. Для классических равновесных моделей т м 0. Для релаксационных моделей, которые учитывают конечное время перераспределения флюидов в капиллярном пространстве, т >  0. В модели [7] в качестве релаксационной переменной используется термодинамический внутренний параметр ф Кинетическое уравнение выбирается таким образом, чтобы гарантировать неотрицательную диссипацию капиллярных сил. Следуя необходимым выводам в [15], управляющие уравнения для релаксационной модели имеют вид

Is + |-b(S,t)=0,(18)

It Ix

b(S,$ =--Р^2, ц = ^,(19)

,^   ki(S)+ ^(S), м Ц1 ’ dS в

S = S + тд = 2S + 7ф - 1,(20)

It=Т [I(1 - S) - «]• граничные и начальные условия:

S (xo,t) = 1; t(xo,t) = 0; S (x,to) = 1; t(x,to) = a,(22)

p где ф - параметр неравновесности, t - безразмерное время, т - безразмерный параметр ре-лаксащш. к - проницаемость с выражениями, приведёнными в Ур. (16). аы|- константы. В данной статье расчёты выполняются с использованием a/Р = 1.

Уравнения численно решаются с использованием TVD-схмы:

' ' = S’ - 1 (B+ 2 - B 2) ’                        - ' „+1 _ СЩУр-S?)                      , Д =     1 + At    , T

B , + 1 = ¹ [b(SfM) + Ь(«У’Ы1)] - ¹»,+1 l^SO - Sf)’            (25)

²      2                                       2     ²

ai+1=max f Ц(sf,t) ’ Ц (8^1’фз+i)) ’(

2         \ Is           Is!

S, = S, - |Ax • AS,,(27)

Sf = S, + 1Ax • AS,,(28)

S t - S.— S_i+1 - S t

AS , = ™od (^~, —д^—) , (29) где верхний индекс п соответствует дискретному времени, нижний индекс j - дискретной координате. Для расчетов в данной статве Аж = 1/200 и At = 1/10000.

4.    Результаты

Графики зависимости насыщенности от координаты S (ж), полученные с помощью сетевой модели, представлены на рис. 4. Для каждой из кривых выполняется условие Qt = const, что означает введение одинакового объёма смачивающей жидкости. Здесь можно заметить, что для более медленных течений амплитуда волны насыщенности больше, тогда как для более быстрых течений амплитуда меньше. Важным преимуществом данной сетевой модели является существование предела для очень высоких скоростей потока Q. Это происходит потому, что приложенная разность давлений значительно превышает капиллярное давление. При медленных течениях мы наблюдаем большую амплитуду и более медленно движущуюся волну, поскольку происходит локальная пропитка областей с низкой проницаемостью. При быстрых течениях капиллярное давление меньше, следовательно, пропитка не происходит; поток в основном проходит через более толстые трубки, которые являются путями с меньшим сопротивлением, и значительное количество несмачивающей жидкости остаётся в областях с низкой проницаемостью; поэтому кривая S (ж) лежит ниже, чем при более медленных скоростях закачки.

Рис. 4. Зависимость иасыщеииости смачивающей жидкости S от безразмерной координаты ж, полученная с использованием сетевой модели: кривая-1 ( Q = 10,t = 100), кривая-2 (Q = 50,t = 20), кривая-3 ( Q = 100, t = 10), кривая-4 (Q = 1000, t = 1) ',Q - скорость потока в безразмерных единицах, t - время после начала закачки в безразмерных единицах

Численное решение с использованием континуальной модели представлено на рис. 5. Все графики построены для одного и того же Q, варьировалось только характерное время релаксации. Для очень быстрой релаксации т = 0.01 (рис. 5а) параметр £ релаксировал очень близко к своему равновесному состоянию 1 — S. В местах, где S = 0, возмущение ещё не достигло этой области, и £ остаётся при своём начальном условии £ = 1. Там же, где S > 0. £ близко соответствует 1 — S.

Однако при т = 1.0 (рис. 56) характеристическое время релаксации становится сравнимым с временем течения. Можно наблюдать, что £ изменился лишь незначительно по сравнению со своим начальным состоянием £ = 1. Это и есть ожидаемое поведение параметра £, который является релаксационным параметром: при изменении условий он медленно адаптируется к новому равновесному значению, причём чем больше т, тем больше времени занимает релаксация.

Рис. 5. Зависимость насыщенности смачивающей жидкости S и параметра неравновесности £ от безразмерной координаты х при моделировании континуальной модели с параметром релаксации для двух различных значений характерного времени релаксации т; различные кривые соответствуют разным моментам безразмерного времени t. Кривая-1 (S(x),t = 0.1), кривая-2 (S (x),t = 0.3), кривая-3 ( S(x),t =0.5), кривая-4 (£(x),t = 0.1), кривая-5 (£(x),t = 0.3), кривая-6 (£(x),t = 0.5)

Рис. 6. Зависимость насыщенности смачивающей жидкости S от безразмерной координаты х по результатам численных вычислений континуальной модели с параметром релаксации, в безразмерный момент времени t = 0.3, для различных значений характериого времени релаксации т: кривая-1 ( т = 0). кривая-2 (т = 0.001). кривая-3 (т = 0.01). кривая-4 (т = 0.1). кривая-5 (т = 1.0). крнвая-G ( т = 10)

На рис. 6 представлены кривые S(x) для различиых значений т, полученные численным решением модели [7]. Кривая-1 соответствует классической равновесной модели при т = 0, которая является решением системы дифференциальных уравнений (14-17). Остальные кривые соответствуют релаксационной модели. Поскольку кривая-2 очень близка к кривой-1, можно утверждать, что при быстрой релаксации решение релаксационной модели стремится к решению Бакли - Леверетта. При более медленной релаксации, когда порядок времени релаксации сравним с порядком времени течения, амплитуда волны насыщенности меньше, а скорость её движения выше. Важно отметить, что кривая-5 и кривая-б почти совпадают, так же как кривая-3 и кривая-4 на рис. 4.

Это означает, что дальнейшее увеличение скорости потока (или времени релаксации) не влияет на распределение насыщенности в потоке. В целом результаты, полученные с помощью сетевой модели (рис. 4), качественно согласуются с результатами континуальной модели [7] (рис. 6).

5.    Заключение

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

Показано, что модель качественно согласуется с классическими результатами и демонстрирует релаксационный характер процессов вытеснения. Установлено, что модель пригодна для моделирования неравновесных явлений в гетерогенных средах и может быть использована для определения параметров континуальных релаксационных моделей типа Кондаурова [6].