CC BY
30
Вытеснение жидкости в пористых автомодельных средах
Г.А. Кузьмин, О.Н. Соболева1
Институт теплофизики СО РАН, Новосибирск, 630090, Россия 1 Институт вычислительной математики и математической геофизики, Новосибирск, 630090, Россия
Численно моделируется вытеснение жидкости из пористой многомасштабной среды. Предполагается, что вытесняющая и вытесняемая жидкости взаимно не растворимы, но обладают схожими физическими свойствами. Для среды принимается статистическая гипотеза масштабного подобия. Проницаемость моделируется логарифмически нормальным распределением. Численно решено уравнение Дарси и найдена скорость фильтрации в кубе на сетке 256x256x256 ячеек. Найдены простейшие статистические характеристики поверхности раздела как функции времени.
1. Введение
Неоднородность пористой среды оказывает существенное влияние на фильтрационные процессы [1, 2]. Естественные пористые среды обычно весьма неоднородны. Крупномасштабные (грубые) детали структуры среды доступны непосредственному наблюдению, например, крупные зоны, на которые можно разбить среду, уверенно выделяемые пласты или прослои. Мелкомасштабные детали в распределении проницаемости и пористости не известны. Их следует учитывать в рамках статистического подхода. Теоретическое определение эффективных параметров требует решения соответствующих задач в средах со случайными полями. В настоящей работе эффективные параметры определяются с использованием усовершенствованных гипотез масштабного подобия Колмогорова для трехмерных сред [35]. Движение жидкости в пористой среде часто изучается с привлечением концепции фракталов [6] и задачи протекания. Применение данных методов требует использования геометрического языка, который недостаточно знаком многим специалистам, работающим с пористыми неоднородными средами. Основные идеи, применяемые в настоящей работе, не выходят за рамки формулировок, содержащихся в цитированных работах Колмогорова. Другим преимуществом применяемого
подхода является то, что дистанция, разделяющая основные параметры и функции теории неоднородных сред от тех, что получаются в экспериментах, представляется более короткой.
2. Масштабное подобие пористой среды
Пористая среда описывается набором полей, которые могут зависеть от пространственных координат и от времени (пористость, проницаемость, коэффициенты в законе Гука и т.д.). Измерение этих полей неизбежно связано с их сглаживанием по пространству и по времени. Зависимость параметров от времени является обычно слабой, поэтому наиболее существенно пространственное сглаживание. Рассмотрим поле проницаемости среды, которое описывает скорость фильтрации жидкости при наличии в ней градиента давления. Пусть несжимаемая жидкость протекает через пористую среду с коэффициентом проницаемости е(х). При малых числах Рейнольдса скорость фильтрации V и давление р связаны законом Дарси V = -е(х)Ур. Условие несжимаемости div V = 0 приводит к уравнению:
Э
дxj
е(х) д- Р(х)
j
= 0.
(1)
© Кузьмин ГА., Соболева О.Н., 2002
Пусть поле проницаемости известно. Это означает, что в каждой точке х мы выполнили ее измерение, прокачивая жидкость через очень малый образец размера 10. Случайная функция пространственных координат е(х) рассматривается как предел проницаемости е(х, 10) (при 10 ^ 0 е(х, 10) ^ е(х)). Зависимость е(х, 10) от масштаба можно рассматривать не как досадную неточность процессов измерения, а как фактор, позволяющий развить новые подходы к исследованию случайно неоднородной среды. Как перейти на более грубую сетку 11 ? Можно было бы просто сгладить полученное поле е(х, 10) по масштабу 11 > 10, но будет ли полученное поле истинной проницаемостью, которая описывает фильтрацию в области масштабов (11, L), где L — наибольший масштаб? Ясно, что это, вообще говоря, будет не так. Чтобы найти проницаемость на более крупной сетке, надо выполнить измерение заново, прокачивая жидкость через более крупные образцы размера 11. Необходимость этой процедуры вызвана тем, что флуктуации проницаемости из интервала масштабов (/0, 11) имеют корреляции с индуцируемыми ими флуктуациями давления. Для их нахождения ниже используется уравнение (1). Поиск закона преобразования эффективной проницаемости при изменении масштаба сетки облегчается в средах, обладающих масштабным подобием. Аналогично Колмогорову [4], Яглому [5], рассмотрим безразмерное поле % равное отношению проницаемостей, сглаженных по двум различным масштабам (I, Г), буквой ~(х) мы обозначаем сглаженную по масштабу I проницаемость е(х, 10): у(х, I, I') = = ~(х, I')/~(х, I). Устремляя Г ^ I, мы получаем поле Ду(х, 1,1у)
ф(х, I) = -
Му
которое определяет все статис-
у=-
тические свойства пористой среды [7]. Полученное соотношение представляет собой дифференциальное уравнение, решение которого дает проницаемость как функцию автомодельного поля. Поле є определяется формулой
є(х, /0) = є0 ехр
■ |ф(х І)у
(2)
Предположим для определенности, что неоднородности более крупного масштаба, чем L, отсутствуют. Мы предполагаем, что корреляционная функция от ф(х, I) изотропна, однородна и масштабно инвариантна,
то есть
Ф(х, у, І, Ґ) = Ф
( (X - у)2 І Л
І2
І
Мы строим простейшую модель, предполагая, что функция ф(х, І) имеет гауссово распределение и дельта-коррелирована по логарифму от масштаба. При х = у мы имеем:
Ф(х, х, І, Г) = (ф(х, І)ф(х, І')) - (ф(х, І))(ф(х, І')) =
= Ф08(1п І - 1п І').
Это предположение соответствует логарифмически нормальной модели для проницаемости. Конкретизацию логнормальной модели мы начнем с замечания, что при любых х, I должно выполняться равенство (е(х, I)^ = е0. Для таких полей, как пористость, это условие естественно следует из их физического смысла. Оно справедливо и для поля проницаемости, так как сглаживание по большим объемам, согласно эргодичес-кой гипотезе, эквивалентно статистическому осреднению. Для гауссова распределения это условие приводит к равенству Ф 0 = 2^ ф.
3. Численное моделирование фильтрации жидкости
Решается задача о фильтрации в неоднородной пористой среде в кубе с ребром L0. На гранях у = 0 и у = L0 задается постоянное давление
р(X у> г)\у=0 = р(X у> г)\у=L = Р2, Р1 > Р2.
Давление на всех прочих гранях куба задается линейной зависимостью р = р1 + у(р2 - рО/^ . Основной фильтрационный поток направлен по оси у. Флуктуации пористости вызывают флуктуации величины и направления скорости фильтрации.
Для численного расчета течения выполним в (1) переход к безразмерным переменным. Все длины измеряются в единицах L0, за единицу разности давления выбирается р1 - р2. Проницаемость измеряется в единицах е0. Таким образом, достаточно решить задачу в единичном кубе с единичным скачком давления при е0 = 1. Вначале выполняется генерация поля проницаемости. Интеграл в формуле (2) заменяется конечно-разностной формулой, в которой удобно перейти к логарифмам по основанию 2:
єІ (х) = єо ехр
1 1о§21 -— {ф(х, т)dт
І0§2 Є
1о§21
= єо2
-Еф(х, Т )8т
(4)
где І = 2Т, 8т — шаг дискретизации по логарифму масштаба. Мы рассматриваем 8т = 1, тг- = -1, -2, ..., -8 и сетку 256x256x256 по пространственным переменным. Используется корреляционная функция
(ф(х,І), ф(У, 0)с =
Ф0
1п2
ехр
(х - у)2
о 2т
8(т-т),
0.2 0.4 0.6 0.8 х
Рис. 1. Изолинии модуля скорости в сечении £ = 0.5
где константа Ф 0 = 2^ф должна быть выбрана из экспериментальных данных для пористых природных сред. Центральные моменты отмечаются значком «с» внизу после закрытия угловой скобки. Согласно работе [2] Ф0 ~ 0.3. Структура корреляционной матрицы позволяет представить ее в виде прямого произведения четырех матриц более низкой размерности и применить для численного моделирования алгоритм «по строкам и столбцам» [8]. Дельта-коррелированность по логарифму масштаба означает, что на каждом масштабе I поле ф(х, I) генерируется независимо. Общий показатель степени в (4) складывается из суммы по статистически независимым слоям. В настоящих расчетах верхние два слоя и нижние три слоя оставлены незаполненными, то есть поле ф на них равно нулю. Пустые два верхние уровня означают, что масштаб самых крупных пульсаций проницаемости L = 1/8. Это позволяет при практическом осреднении приближенно заменить вероятностные средние усреднением по пространству. То, что самые мелкие три уровня оставлены пустыми, обусловлено требованием, чтобы наша разностная задача хорошо аппроксимировала бы уравнение (1) на всех масштабах. Для приближенного расчета используется ограниченное количество слоев. У нас число этих слоев равно 4.
Затем численно решается сеточный аналог безразмерного уравнения (1). Используется итерационный метод в сочетании с быстрым преобразованием Фурье и методом прогонки второго порядка точности [9]. Изолинии скорости приведены для случая, когда заполнены четыре слоя по масштабу с номерами от -6 до -3 (рис. 1). Далее мы вычисляем суммарный расход, используя прямое численное решение уравнения (1). Дан-
ная величина получена при заполненных четырех слоях по масштабу, затем пульсации самого мелкого масштаба обнуляются, и вычисляется расход. Затем оставляются свободными два самых мелкомасштабных слоя и повторяются вычисления расхода и так далее. Результат приведен на рис. 2 в двойном логарифмическом масштабе. График демонстрирует степенную зависимость логарифма расхода
от логарифма минимального масштаба т. Наклон ~ 0.1 (линия 1) к оси абсцисс при Ф0 ~ 0.3 и наклон 0.2 (линия 2) при Ф0 ~ 0.6. Хотя мелкомасштабные пульсации выбираются так, что они не влияют на среднюю по пространству величину 8, полный расход через среду зависит от того, насколько полно учитываются мелкомасштабные пульсации. Данный результат показывает необходимость учета мелкомасштабных пульсаций при расчете крупномасштабной компоненты скорости. Приемлемая формула учета мелкомасштабных пульсаций для логнормальной проницаемости приведена в [10, 11].
4. Вытеснение жидкости в пористой среде
Выполненные расчеты могут послужить для изучения вытеснения одной жидкости с помощью сходной с ней, но другой жидкости. Поверхность раздела помечается пассивными частицами, которые в начальный момент времени занимают некоторые заданные положения, затем перемещаются стационарным полем скорости. Поскольку жидкости физически одинаковы, то их скорость фильтрации удовлетворяет уравнению (1), которое было численно решено в предыдущем параграфе. Движение точек поверхности раздела описывается уравнением [1]
Нх I
«(х)— = Чх), х (=0 = х0; , i = 1 ..., Ы, (6)
где т(х) — пористость; V — скорость фильтрации, приведенная на рис. 1. Случайное поле пористости коррелирует с проницаемостью. Пористость моделируется
а
1
Рис. 2. Зависимость О (логарифма полного расхода через среду) от минимального масштаба т
Ур
Р1 - р 2
0 0 0 0 Рис. 3. Поверхность раздела двух жидкостей в момент времени Ц (а) и в момент времени 12 (б), 12 > t1
б
1
таким же образом, как это было сделано выше для проницаемости
т(х, /0) = т0 ехр
-|| Х(х, /) |у
(7)
где вместо ф(х, /) (см. (2)) используется другая случайная безразмерная функция %(х, /). Появление абсолютного значения | х(х, /) | мотивируется дополнительным ограничением на пористость т(х, /0) < 1, которое следует из ее определения. Если бы знак абсолютной величины не использовался, отрицательные значения флуктуирующего поля привели бы к нарушению ограничения на т(х). Чтобы обеспечить справедливость выполнения условия (т(х)) = т0, мы вычисляем среднее значение интеграла (7) для одного слоя в точке х, которое используем как нормировочный множитель при вычислениях с несколькими слоями. Статистическая связь пористости и проницаемости фиксируется корреляци-
онной функцией ^ф(х, /)х(у, /Эволюция поверхности раздела во времени приведена на рис. 3.
Рис. 3 показывает, что форма поверхности раздела сильно усложняется во времени. Точки поверхности движутся подобно броуновским частицам. Как показано на рис. 4 (в координатах у-), среднее движение поверхности раздела происходит с некоторой постоянной скоростью.
На рис. 5 приведен средний квадрат отклонения точек поверхности от среднего уровня 0. Используются двойные логарифмические координаты. Если отбросить входной и выходной участки, где ощущается влияние граничных условий, увеличение квадрата толщины происходит по степенному закону с показателем ~ 1.73. Значения параметров: Ф0 ~ 0.2 (для проницаемости), Е 0 ~ 0.3 (для пористости), где Е0 — константа, аналогичная Ф0 в формуле (3), при умеренных значениях взаимной корреляции. Величина 1.73 превосходит классическое броуновское значение 1.0, что иллюстрирует
Рис. 4. Зависимость {у(£)^
Рис. 5. Зависимость 0 от 1пt, где 0 = 1п^(у -(у^)2^
существенное различие рассматриваемой диффузии поверхности с классическим случаем.
5. Заключение
Простой взгляд на приведенные графики обнаруживает их сходство с результатами численных экспериментов с использованием методологии мультифракталов. В рамках подхода, изложенного в настоящей работе, получатся мультифракталы, если минимальный масштаб І0 устремить к нулю. В настоящей работе минимальный масштаб остается конечным, поэтому сингулярности отсутствуют. Канторовы множества не возникают, и весь анализ не выходит за пределы аппарата дифференциальных уравнений и теории случайных процессов. Основными объектами теории являются поля, свойства которых могут быть (хотя бы в принципе) напрямую измерены экспериментально.
Конвективная диффузия описывается симметричным тензором второго порядка Бу (І). Тензор Бу (І) необходим при разработке феноменологических и под-сеточных моделей вытеснения жидкости. Мы выполнили вычисление всех компонент Бу (І) и выяснили, что поперечная диффузия мала по сравнению с продольной (по отношению к направлению течения). Этот результат влияет на выбор моделей распространения примеси и вытеснения жидкости. Мы планируем рассмотреть его подробнее в отдельной работе.
Простейшей моделью, которая рассмотрена в рамках настоящего подхода, является масштабно инвариантная логарифмически нормальная модель несущей среды. Логарифмически нормальная модель критиковалась по нескольким причинам. Например, логарифмически нормальное распределение не удовлетворяет условиям теоремы Карлемана, из-за чего оно не определяется своими моментами. Имеются и работы, в которых приводятся аргументы, что при правильном применении логнормальной модели больших проблем не возникает [12]. Гипотеза подобия Колмогорова, сама по себе, не предполагает логарифмической нормальности. Мы использовали логнормальную гипотезу для просто-
ты и потому, что это предположение не ведет к каким-либо осложнениям в нашем случае. Другая очевидная возможность — изучить другие устойчивые распределения [13, 14].
Главными результатами работы следует считать полученные степенные законы фильтрации и диффузии в пористой автомодельной среде.
Работа выполнялась при поддержке РФФИ.
Литература
1. Швидлер M-И. Статистическая гидродинамика пористых сред. -М.: Недра, 1985. - 287 с.
2. Sahimi M. Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata, and simulated annealing // Reviews of Modern Physics. - 1993. - V. б5. - P. 1393-1534.
3. КолмогоровA.H. О логарифмически нормальном законе распреде-
ления размеров частиц при дроблении // Докл. АН СССР. - 1941. -Т. 31. - С. 99-1G1.
4. Kolmogorov A.N. A refinement of previous hypotheses concerning the
local structure of turbulence in a viscous incompressible fluid at high Reynolds number // J. Fluid Mech. - 19б2. - V. 13. - P. 82-85.
5. Moнuн A-С., Яглом A.M. Статистическая гидромеханика. - М.: Наука. - 19б7. - Т. 2. - 72G с.
6. Mandelbrot B. The fractal geometry of nature. - New York: Freeman and Company, 1983. - 45G p.
7. Kuz’min G., Soboleva O. Conformal symmetric model of the porous media // Applied Mathematics Letters. - 2GG1. - V. 14. - P. 783-788.
8. Ogorodnikov V.A., Prigarin S.M. Numerical modeling of random processes and fields: algorithms and applications. - Utrecht-Tokyo: VSP, 199б. - 245 p.
9. Mapчук Г.И. Методы вычислительной математики. - М.: Наука, 1989.
1G. Kuz’min G.A., Soboleva O.N. The renormgroup model for fluid flow through the fractal porous media // Extended Abstracts of Reservoir Rock, EAGE/SEG Res. Workshop, Pau, 2GG1.
11. Кузьмин r.A., Соболева O.H. Подсеточное моделирование фильтрации в пористых автомодельных средах // ПМТФ. - Т. 43. - №4.-С. 115-12б.
12. Moлчaн r.M. Турбулентные каскады: мультифрактальные характеристики // Вычислительная сейсмология. - 1997. - Т. 29. - С. 155-1б7.
13. Shertzer D., LoveJoy S. Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes // J. Geophys Res. - 1987. - V. 92. - No. D8. - P. 9б9З-9714.
14. Bouffadel M.C., Lu S. et al. Multifractal scaling of the intrinsic permeability // Water Resource Research. - 2GGG. - V. Зб. - No. 11. -P. 3211-3222.
