Моделирование двумерной фильтрации меченой жидкости методом трубок тока Текст научной статьи по специальности «Математика»

Лялин В.Е.; Сидельников К.А. МОДЕЛИРОВАНИЕ ДВУМЕРНОЙ ФИЛЬТРАЦИИ МЕЧЕНОЙ ЖИДКОСТИ МЕТОДОМ ТРУБОК ТОКА

Постановка проблемы

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

Уравнения многофазной многокомпонентной фильтрации

Движение меченой жидкости через пористую среду можно рассматривать как частный случай многофазной многокомпонентной фильтрации флюидов, поэтому для начала запишем общие уравнения сохранения для изотермического потока i-го компонента [0; 0]:

Я пр пр пр

+ v • X - v • =о, /=î.. .и„, 11 )

ol j=1 j=i j=i

где ф - пористость (объем порового пространства в единичном объеме пористой среды); - мо-

лярная доля i-го компонента в j -й фазе; р^ - молярная плотность фазы j ; Sj - насыщенность j -й фазой (объем, занимаемый j -й фазой, в единичном объеме порового пространства); й- - вектор скорости фильтрации j -й фазы в случае многомерного потока; - тензор гидродинамической дисперсии I -

го компонента в j -й фазе; п - общее число фаз; пс - общее число компонентов в системе.

Кроме того, для уравнения (1) устанавливаются дополнительные ограничения [0; 0]:

пР пс

j=l-«P- (2)

j=1 /=1

Воспользуемся обобщенным уравнением Дарси (как частного случая закона сохранения импульса) предложенного Маскетом для моделирования скорости фильтрации [0; 0]:

(VPJ + = l-”Pr ( 3 )

О

где k - тензор абсолютной проницаемости (в большинстве практических задач можно (или необходимо) предположить, что он диагональный); krj - относительная проницаемость для j -Й фазы; //. - вязкость j -й фазы; Р. - давление в j -й фазе; р . - массовая плотность j -й фазы; g - вектор уско-

рение свободного падения.

Нижний индекс для давления в уравнении (3) означает, что в разных фазах давления имеют разную величину, другими словами фазы разделены кривой поверхностью раздела с ненулевым межфазным натяжением. Связь между этими давлениями всегда представляется через функцию капиллярного давления.

Строго говоря, под гидродинамической дисперсией понимается параллельное существование молекулярной диффузии и механической дисперсии, которые являются независимыми физическими явлениями. Тем не менее, математические представления этих двух механизмов очень похожи, а результаты их влияния обычно эквивалентны друг другу [0]. Тензор гидродинамической дисперсии D^ i-го компонента в j -й фазе включает как диффузионные так и дисперсионные члены:

D^D.+D,, (4)

где D(( - тензор молекулярной диффузии i -го компонента в j -й фазе; D^ - тензор механической дисперсии

i-го компонента в j -й фазе. Типичное выражение для компонентов тензора механической дисперсии D^ имеет следующий вид [0; 0]:

/г, \ (а> lUJ s (а‘ )j ~(а‘ j W Ч

( mn \j ф8^ mn + ф8^ (lUm )j (Un )j , (5)

где материальные константы (^t)j и (^l)j - дисперсности среды в поперечном и продольном направлениях потока j -й фазы; 8тп - дельта-функция Кронекера ( ^Ш1= 1 при m = п и 8шп =0 в остальных случаях); m, п = 1, 2, 3 .

Уравнения конвективной диффузии

Выведем уравнения, описывающие фильтрацию меченой жидкости, из общих уравнений для многофазной многокомпонентной фильтрации записанных выше. Для этого предположим, что существует однофазный поток жидкости (и =1, iS^ = 1 , Pi=P, Рс— 0, u\=u, kri= 1, D^=Dj)f которая является несжимаемой с

ПОСТОЯННЫМИ ПЛОТНОСТЬЮ ( р\ = Р , рт\ — Рт ) и вязкостью ( = ¡л ) . Также будем считать, что пласт, че-

рез который протекает жидкость, несжимаем (ф = const ); отсутствуют гравитационные силы ( g = 0); в системе присутствуют только два компонента ( пс =2 , ®n=®i ), причем первый из них является индикатором. В результате уравнения (1), (2) и (3) упростятся и примут вид

^+v(®,.«)-^v(d;-v®,.) = 0, / = 1,2, (6)

Т,ті =1' (7)

/=1

й =----УР . (8)

м

Мы получили систему из четырех уравнений с четырьмя неизвестными: щ , и и Р • Тем не ме-

нее, систему уравнений (6)-(8) можно записать в более удобной для решения форме. Если просуммировать уравнение (б) по всем компонентам, то, учитывая равенство (7), получим Уй = 0 . (9)

Поскольку молярную долю 2-го компонента легко вычислить на основании равенства (7), то система из четырех уравнений (6)-(8) преобразуется в систему всего из двух уравнений с двумя неизвестными СО и Р :

ф^-+й-Ут-ф-У(Т)'-Ут) = 0, (10)

-Ур| = 0 , (11)

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

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

(ЯуIРы , тогда объемная доля і -го компонента в у -й фазе определяется выражением

С, = У^ . (12)

1 = 1

Молярная плотность фазы может быть записана в виде

( пс Л

Ра

-1

у/. с | (13)

Сравнение уравнений (12) и (13) дает РсС а =Р№ • (14)

Подставив уравнение (14) в уравнение (1), и проделав те же упрощения и действия, что и при выводе уравнения (10), получим уравнение конвективной диффузии, но уже для объемной доли компонента индикатора: с)с

ф— + й-Ус-ф-У(1У-Ус) = 0 . (15)

Здесь так же, как и для уравнения (10), будет выполняться равенство (9), поскольку ^ с. = 1 .

^^! = \ 1

Следовательно, можно записать и решать обобщенное уравнение конвективной диффузии:

ЯС

ф-^ + й • УС - ^ • У(Б' • УС) = 0 , (16)

где неизвестная переменная С является долей (молярной или объемной) индикатора в жидкой фазе. Выражение для коэффициента гидродинамической диффузии примет вид (см. уравнение (4)

1У = 1) + 1) , (17)

где компоненты коэффициента механической дисперсии индикатора будут определяться на основании уравнения (5) следующим образом: а*и , а, -а*

Пшп =— 5шп +—-------ишМп • (18)

ф фи

Метод трубок тока

Метод трубок тока основан на двух ключевых идеях: (а) генерация линий/трубок тока для некоторой

интересующей области и (Ь) отображение одномерного решения, полученного для каждой трубки тока. Поскольку для несжимаемого потока выполняются равенства (9) и (11) , линии/трубки тока могут быть найдены путем решения уравнения для функции тока ^ [0]:

д_

дл

к дл

V

д 1 дТ . Л

-----1---------| = 0 . (19)

ду I кх ду

Ключевая идея использования трубок тока для моделирования двумерной процессов вытеснения состоит в рассмотрении каждой трубки тока как одномерной системы. Хиггинс и Лейтон (см. [0]) показали, что для отображения одномерного решения, полученного для трубок тока, оно должно масштабироваться по объему. Представление каждой трубки тока как одномерной системы автоматически свяжет с ней определенный поровый объем - доли общего порового объема пласта. По определению, трубка тока описывает объемный расход. Таким образом, для каждой трубки тока можно использовать обобщенное представление безразмерного времени [0]:

ів = (2°)

ур

где q - расход через трубку тока, равный разности Д^=^в — , где нижние индексы А и В обо-

значают граничные линии тока; Ур - некоторый произвольный поровый объем, используемый при масштабировании.

Аналогично к каждой трубки тока может быть привязана безразмерная длина [0]:

.у

ф{с)йс

" , (21)

" Ур

где А - площадь трубки тока как функции координаты С вдоль трубки тока. «Наилучшим» выбором для величины Ур ясно является [0] тТ У?

У =1Т' (22)

1\зт

где У^т - общий объем пор системы.

Решение одномерного уравнения конвективной диффузии

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

д*в _ Ч

~дГ = У? (23)

и

дхп фА {у)

—^ = _> ’ , (24)

ду Ур

где у - координата при движении вдоль оси трубки тока.

Из равенства (23) легко получить временное преобразование:

дс__ч_ дс

Зг УР ЗГд

Выводы координатного преобразования для вектора скорости филвтрации й и тензора дисперсии О’ индикатора внутри жидкости существенно различаются. Из определения градиента функции известно, что ^ дС -УС = --е„, (26)

где е - единичный вектор в направлении $ . При движении вдоль трубки тока векторы й и е будут коллинеарными, и координатное преобразование для вектора скорости будет иметь вид (см. равенство

(24)

_ ^ _ (дС А /- ЛА дС

и • УС = и ---е„ = (и • е.)---= иф—-----. (27 )

у К дз УР дхв

Для простоты предположим, что влияние молекулярной диффузии незначительно и значением коэффициента 1)тп можно пренебречь, тогда 1Утп — Птп и для двумерного пространства на основании уравнения (18) можно записать равенство

В' = Л

Лі 2 а/

—-------------------‘-их —----------------и и у

-а, 2 ^иу

(28)

Исполвзуя уравнение (26), получим Т>' ■УС = (Т>'-ё\— = б— , (29)

где результат скалярного произведения тензора (28) на вектор е справа:

3 = ^-й . (30)

ф

Из этого следует, что трубка тока как одномерная система не может учесть тензорную природу О или (когда действием молекулярной диффузии можно пренебречь). Поскольку продольная дисперсность

а является величиной постоянной, то на основании уравнения неразрывности (9) и уравнения (29), будет справедлива запись

у(т>' ■ус) = у[ь—1 = £>.уГ—1 = .

v ’ У &) I. а*) &2 &2

(31)

Наконец используя уравнение (24) легко вывести координатное преобразование для тензора дисперсии:

тлд2С ША д Г А дС) ^2 А д ( лдС )

Б—г- = Бф=------ф---------1 = Бф2^----------1 А--I . (32)

Эу УР ддв ^ У дхв 1 Ур дхв ^ схв )

Подстановка уравнений (25), (27) и (32) в уравнение (16) дает

q дС А дС 3 А д І дС

Ф=----------------+ и^—-------------------------------------------------Вф ^-1 А- = 0 . (33)

Ур діВ Ур длВ Ур длВ | длВ ^

а -а

аи а

Известно, что для каждой трубки тока выполняется д=и•А , тогда после несложных преобразований уравнение (33) примет вид

дС дС ф2в д ( лдС . „

+-----?-=---------------1 А-1 = 0 . (34)

дґв дхв иУР дхв I дхв

В общем случае величина площади поперечного сечения трубки тока А зависит от координаты 5 , а значит и от безразмерной длины трубки х^ . Однако предположим, что влияниям величины этой производной на результат решения уравнения (34) можно пренебречь:

£¿4 Л „

---= 0, (35)

с1хв

тогда уравнение (34) упроститься:

дС дС ф2БА д2С „

— +--------2-=------ = 0 . (36)

дгв дхв иУр дхв

Для пласта с площадью поперечного сечения Ар и длиной Ь общий поровый объем Ур = фАрЬ , откуда на основании уравнения (22) запишем

тТ фА„Ь

Ур =. (37)

Подставив уравнение (37) в уравнение (36), получим

дС дС фв КЗТА д2С

дів дхв иЬ Ар дхв 2

=0. (38)

Если число трубок тока в системе достаточно велико и длина пласта Ь существенно превышает его толщину Н , то вполне достоверно предположить, что

А\ТА Л 1

5Т =1 . (39)

Лг

В результате уравнение (38) примет вид

дС дС фв д2С „

----+---------------- = 0 . (40)

д1в дхБ иЬ дхв 2

Введем обозначение для так называемого безразмерного числа Пекле п иЬ

Ре =---, (41)

фВ

представляющего собой отношение характеристического времени дисперсии (Ь2/В ) к характеристическому времени конвекции фЬ/и . Подставив уравнение (30) в уравнение (41), получим еще одну формулу для числа Пекле: г, Ь

Ре = — . (42)

а,

В результате уравнение (40) окончательно запишется в следующем виде:

дС дС 1 д2С „

— +---------------- = 0 . (43)

дгв дхв Ре дхв

Легко заметить, что при большом значении числа Пекле влияние дисперсии мало и распространение

индикатора внутри жидкости осуществляется главным образом за счет конвекции, тогда уравнение (43)

примет вид:

дС дС „

----+-------= 0 . (44)

дГв дхв

Решение уравнения (44) с условиями соответствующими задаче Римана представляет собой простое

движение волны с единичной скоростью:

["1, для хв < 1^в

(45)

0, для хв > хв

Ясно, что уравнение (44) является предельным случаем, т. к. его решение (45) представляет собой резкий фронт вытеснения (поршневое вытеснение), движущийся с единичной скоростью. Тем не менее, уравнение (44) остается важным, потому что его решение является эталонным, с которым сравниваются другие численные решения.

Если число Пекле постоянное (не зависит от состава и характеристик пористой среды), то уравнение (43) может быть легко решено с использованием преобразования Лапласа. Для условий задачи Римана решением является [0; 0]

( + \\

• (46)

с (хв, в )=іег&

Ограничения координатного преобразования.

Для преобразования уравнения (16) в его одномерное представление, решение которого затем отображается на каждую трубку тока, были сделаны два основных допущения касательно коэффициента гидродинамической дисперсии (см. формулу (30) и геометрии трубок тока (см. выражения (35) и (39).

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

что а ¡и »В (см. также [0] ) .

В своей работе Тил [0] привел ссылки на исследования, в которых показано, что с увеличением скорости фильтрации начинает преобладать диффузия в продольном направлении, при этом величины двух коэффициентов диффузии могут различаться на несколько порядок. Этот вывод справедлив также и для фильтрации в масштабе месторождения. В большинстве случаев это подтверждается промысловыми данными). Таким образом, перенося решения на рис. 3 в масштаб месторождения, оно будет не сильно отличаться от решений, где бы учитывалась диффузия в поперечном направлении.

При выводе уравнения (43) были сделаны два предположения относительно геометрии трубок тока, связанные с площадью поперечного сечения трубки тока а (*) как функции координаты * вдоль трубки

тока и производной этой функции относительно безразмерной длины трубки тока хв(*) .

В общем случае функция А (*) может принимать разнообразные значения, причем графики этих функций

для разных трубок тока существенно различаются, а значит отношение (39) далеко от единичного значе-

ния. Разброс значений Ы5ТА/Ак , где А - средняя площадь поперечного сечения трубки тока, относительно единицы также достаточно большой и в некоторых случаях может достигать двукратной относительно единицы величины. Введем модифицированное число Пекле ' —\-1

Ре = Ре .\Ы^ГА

I Ая

которое будет иметь какое-то свое значение для разных трубок тока. В результате вместо решения (4 6) одинакового для всех трубок тока, необходимо использовать различные одномерные решения для отображения на каждую трубку тока

( Л

с(хв, 1в) = ^егй; -—; ^ в' . (47)

Ре ( хв (в )

2^/(B"

Легко проверить, что даже при двукратном отличии модифицированного числа Пекле Ре от первоначального значения Ре , максимальное значение разности функций (46) и (47) по абсолютной величине не превышает 0,084 и при отображении практически не заметно. Значит, подобная поправка не меняет существенно вид решения и предположение (39) вполне допустимо.

Запишем производную СА( 5)/Зх^ в следующем виде:

сіА( 5) Ур 1 СА( 5)

СX ф А( 5) І5

(48)

Влияние производной (48) на результат решения одномерного уравнения конвективной диффузии выражается в появлении дополнительного конвективного члена (см. уравнение (34)

1 ф2в СА(5)дС

иУр (Схв )дхв

(49)

Другими словами, случайные изменения площади поперечного сечения приводит к колебанию безразмерной скорости конвекции относительно единичного значения. На базе формул (41), (42) и (48), уравне-

ние (49) примет вид

1 ЛА(*)^ дС

1 -а, • ———— —. (50)

А(*) Л* ^дхв

Известно (см. [0]), что щхЬ , где Ь - средняя длина порового канала, и коэффициент пропорциональности зависит от дисперсии средней проводимости среды в различных направлениях. В частности асимптотическое выражение дисперсности пористой среды в продольном направлении для поля проницаемости

с малой дисперсией а= 0,5^^кЬ [0] . Однако для того, чтобы модель трубки тока могла трактоваться как система, удовлетворяющая ограничениям законов Фика, должно выполняться неравенство L«ylc«L, (51)

где Хс - корреляционная длина поля проницаемости, Ь - размер области в основном направлении потока. Если порядок величины Хс обычно 10 1 (см. [0]), то на основании левого неравенства уравнения (51) порядок Ь не должен превышать 10 2 . С другой стороны порядок выражения

1 СА( 5 )

А( 5 ) І5

обычно не превышает 10 1 для «мгновенных» значений и 10 3 для средних. В результате влиянием дополнительного конвективного члена вносимого ненулевым значением производной (48) ввиду его малости можно пренебречь.

Учебный пример

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

то доля индикатора сразу определяется по формуле (45) или (46) в заданный момент времени для

всех точек отображения. Отображение решения (46) на каждую трубку тока немного отличается от случая чистой конвекции. Для нахождения решения соответствующего определенному моменту времени (в необходимо явно вычислять выражение в правой части формулы (4 6) для этого момента. Пример двумерных

решений для неоднородной области для момента времени (в = 0,55 и трех различных значений числа Пекле показан на рис. 3.

Точность двумерного решения будет зависеть от числа трубок тока используемых для отображения одномерного решения, а также от числа точек отображения для перемещения решения от трубок тока обратно на обычную прямоугольную сетку. Ясно, что чем больше трубок тока, точек отображения и точек, образующих граничные линии тока, тем точнее решение. С другой стороны должен существовать предел необходимого количества трубок тока для представления двумерного решения. Если каждая трубка тока будет иметь максимальную ширину, которая меньше, чем размер блока основной прямоугольной сетки, то положение фронта из-за большого числа трубок тока будет потеряно для операции усреднения на уровне сеточного блока.

Рис. 1. Линии ницаемости

(трубки) тока и линии равного давления для недородного пласта с анизотропией про-

Рис. 2. Пример решения без учета диффузии в момент времени ^ = 0,55

Ре = 10-

Ре = 103

Ре = 10

Рис. 3. Решения уравнения (16) для трех чисел в момент времени = 0,55

ЛИТЕРАТУРА

1. Bear, J., Dynamics of fluids in porous media, Dover, New York, 2005. Reprint. Originally published: American Elsevier, New York, 1972.

2. Juanes, R., Displacement theory and multiscale numerical modeling of three-phase flow in po-

rous media, Ph. D. Thesis, University of California, Berkeley, California, 2003.

3. Orr, F. M. Jr., Theory of gas injection processes, Stanford University, Stanford, 2005.

4. Thiele, M. R., Modeling multiphase flow in heterogeneous media using streamtubes, Ph. D. The-

sis, Stanford University, Stanford, California, 1994.