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

ДОБЫЧА

УДК 622.276

методика измерения кривых капиллярного давления дренирования по цифровым микротомографическим изображениям пористых сред

Т.Р. Закиров

к.ф.-м.н., научный сотрудник timurz0@mail.ru

А.А. Галеев

к.ф.-м.н., доцент akhmet.galeev@kpfu.ru

ФГАОУВО Казанский (Приволжский) федеральный университет, Институт геологии и нефтегазовых технологий, Казань, Россия

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

материалы и методы

Для математического моделирования двухфазных течений в поровых каналах используются решеточные уравнения Больцмана. Явления на границе раздела двух фаз описываются при помощи градиента цифрового поля. В качестве образцов исследования используются цифровые модели типичных при тестировании песчаников Береа и 1_\/60А.

Ключевые слова

математическое моделирование, кривые капиллярного давления, рентгеновская микротомография, пористая среда, решеточные уравнения Больцмана

Фильтрационные свойства пористых сред, такие как коэффициент абсолютной проницаемости, кривые капиллярного давления и относительные фазовые проницаемости, являются одними из немногих характеристик горных пород, которые нельзя оценить по геофизическим данным. Их измерение, как правило, проводится экспериментально в лабораторных условиях и предполагает ряд трудоемких этапов, среди которых подготовка флюидов и поверхности поровых каналов, а также длительное время проведения экспериментов на слабопроницаемых образцах. В качестве альтернативного подхода в последние годы интенсивно разрабатываются методики проведения численного эксперимента на основе цифрового микротомографического изображения керна (computational rock physics). Если измерение коэффициентов абсолютной проницаемости при помощи различных моделей (уравнения Навье - Стокса, решеточные уравнения Больцмана) является уже хорошо исследованным и зарекомендовавшим себя методом [1, 2], то измерение кривых капиллярного давления и относительных фазовых про-ницаемостей содержит ряд трудностей, среди которых: математическая постановка задачи и длительное время расчетов. На сегодняшний день известен подход, согласно которому кривые капиллярного давления измеряются по модели поровых сетей или pore network models (PNM) [3], в которых исходному цифровому образу порового пространства ставится в соответствие модель из сферических пор и соединяющих их цилиндрических капилляров. По такой модели с заданной геометрией пор и капилляров проводится вычислительный эксперимент, и капиллярное давление оценивается по формуле Лапласа. Все более широкая доступность вычислительной техники с высоким быстродействием и усовершенствование компиляторов позволяют использовать подход, основанный на проведении вычислительного

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

Постановка вычислительного эксперимента для определения кривых капиллярного давления при дренировании

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

Измерение насыщенности образца смачивающей фазой на каждом режиме дренирования проводится при условии прекращения вытеснения смачивающей и несмачивающей фаз соответственно, а созданный при этом перепад давления, следовательно, является капиллярным.

60

Рис. 1 — Бинаризованная цифровая модель песчаника Береа размером 200x200x200 вокселей, разрешение 3,2 мкм: а - трехмерная проекция; б - срез в плоскости XY, серое - скелет, синее - поровое пространство

Fig. 1 — Binary digital model of Berea sandstone, size 200x200x200 voxels, resolution 3.2 цш: a - 3D view; b - XY plane cross-section view; grey - pore matrix, blue - pore volume

Экспозиция НЕфть газ ноябрь 7 (60) 2017

Плотности смачивающей и несмачиваю-щей фаз во всех вычислительных экспериментах равны 1000 кг/м3, а вязкости —2 мПа*с и 1 мПа*с соответственно. Краевой угол смачивания на границе раздела смачивающая жидкость / несмачивающая жидкость / твердая фаза для данных песчаников, согласно работам [4, 5], равен 0°. Используемая модель предполагает постоянство температуры, химического состава, отсутствие деформации скелета и плотности жидкостей в используемом диапазоне давлений.

математическая постановка задачи. Решеточные уравнения Больцмана

Для математического описания двухфазного течения жидкостей в данной работе используются решеточные уравнения Больцмана [6, 7] (Lattice Boltzmann method, LBM), а явления на границе раздела фаз описываются при помощи модели градиента цветового поля [6]. Выбор решеточных уравнений Больцмана для исследования многофазных течений в цифровом керне обусловлен относительной простотой их численной реализации, а также хорошей адаптируемостью к параллельным вычислительным технологиям, что особенно важно при моделировании на сетках, содержащих, как в нашем случае, несколько миллионов узлов. Данная математическая модель описана во многих работах, поэтому мы не будем давать ее подробное математическое описание, а остановимся лишь на основных базовых положениях.

В рамках LBM течение среды рассматривается с точки зрения динамики ансамбля частиц с заданным конечным числом возможных скоростей. Область течения разбивается сеткой с ячейками, как правило, квадратной или кубической формы. Совокупность данных ячеек составляет решетку (lattice). За шаг по времени At частицы без взаимодействия друг с другом могут совершить один акт перехода между соседними узлами решетки. В качестве переменных, описывающих состояние системы в каждом узле сетки, используются одночастичные функции распределения f (r, u, t). Данная функция показывает долю частиц в момент времени t, находящихся в окрестности точки r (x, y, z) с координатами от x до x+Ax, от y до y+Ay, от z до z+Az и со скоростями в диапазоне от u (u, u, u) до u (u+Au, u+Au, u+Au).

Для двухмерной области течения используется дискретный набор скоростей D2Q9, который задается следующим образом: et = с(0,0); e2 = с-(1,0); e3 = c-(0,1); et = c-(-1,0); e5 = c-(0, -1); e6 = c-(1,1); e7 = c-(-1,1); e8 = c-(-1, -1); e9 = c-(1, -1), где — At / At — базовая скорость (л! — шаг сетки).

Каждому вектору скорости из заданной совокупности ei (/=1,9) ставится в соответствие функция распределения f (r, t), зависящая только от t и r.

Описание динамики ансамбля частиц каждого из флюидов проводится в несколько этапов. Первый этап является этапом переноса (streaming step). На данном этапе частицы за время At смещаются в соседние узлы, причем направление скорости движения не меняется. На втором этапе рассматривается процесс столкновения частиц (collision step), в результате которого функция распределения частиц стремится к равновесному состоянию. На третьем этапе описывается взаимодействие флюидов друг с другом на границе раздела, а также с твердой фазой. Эволюция функции распределения каждого флюида во времени и пространстве описывается при помощи уравнения (1):

f,l(r + e,M,t + M) = f,k(r, )( + (Of (г, i))1 + (21) (f, i))2 (1)

где k = 1,2указывает на вид флюида, т.е. на смачивающую и несмачивающую фазы соответственно.

В зависимости от формы представления оператора столкновения (Q f)' , различают SRT (single relaxation time) модель [6] и MRT (multitriple relaxation time) модель [9]. Авторами работы [9] показано, что данные, полученные с использованием MRT-модели, имеют более точное совпадение при сравнении с результатами аналитического решения на примере течения Куэтта. По данной причине в работе используется MRT-модель.

Параметр релаксации f является определяющим в SRT и MRT моделях. Он служит для управления кинематической вязкостью pk флюида и связан с ней соотношением (2):

M

2тк -1 j А/2 6 I At

(2)

Уравнения в методе 1_ВМ решаются в переменных «плотность-скорость». Макроскопические плотность и компоненты скорости каждого флюида в ячейках вычисляются по формулам (3) и (4):

Рис. 2 — Бинаризованная цифровая модель песчаника LV60A размером 200x200x200 вокселей,разрешение 10,0 мкм: а - трехмерная проекция; б - срез в плоскости ХУ;серое - скелет, синее - поровое пространство

Fig. 2 — Binary digital model of LV60A sandstone, size 200x200x200 voxels, resolution 10.0 цш: a - 3D view; b - XYplane cross-section view; grey - pore matrix, blue - pore volume

pk (r,t)=t fk (r, t) (3)

uk (r,t) = p t e f (r, t) , (4)

Давление pk в методе LBM, создаваемое каждым флюидом, связано с его плотностью следующим соотношением: pf [6, 7].

Для описания явлений, прои3сходящих на границе раздела флюидов, в данной работе используется Color-field метод [8], состоящий из нескольких этапов:

1) Расчет градиента цветового поля g, компоненты которого вычисляются по формуле:

g(r,t) = fiel(fl2(r + el&t,t)-fl1(r + el&t,t)). (5)

¿-1

Традиционно одному из флюидов ставится в соответствие красный цвет, а второму — синий.

2) Описание эффектов поверхностного натяжения на границе раздела флюидов:

( "f )2 = A g (2 • cos2 (а) -1), (6)

где A — параметр, управляющий поверхностным натяжением, at—угол между вектором g и направлением ei.

3) «Recoloring» step — модификация функции ff после решения уравнения (1):

С/)2)" =—/,+/

Р

Р Р

fr

• соs(а, ) ,

• cos (•••• S >

, f.sc — равно-

p p

где р = рг + P , f = f. + f. весная функция распределения [6, 7], рассчитанная при плотности р и скорости, равной нулю. Параметр ß управляет толщиной границы раздела флюидов. В данной работе его значение составляет 0,8 (и не может превышать единицу). Уменьшение данного значения приводит к «размазыванию» межфазного фронта, а увеличение — к нарушению изотропности пространства в задаче о формировании сферической капли (тест Лапласа) [10].

В качестве граничных условий на твердых и внешних непроницаемых границах области течения используются условия «bounceback» [6, 7], являющиеся аналогами условий непротекания и прилипания жидкости в классической постановке задачи. На входной и выходной границах считаются известными давление жидкости и нормальные к границе компоненты скорости, которые равны нулю. В LBM такие условия задаются при помощи соотношений Зоу и Хе [11].

Модель была реализована в среде разработки Visual Studio 2012 на языке Intel Fortran. Расчеты производились на вычислительной станции, содержащей два процессора Intel Xeon, каждый из которых имеет по 20 логических ядер. Для ускорения расчетов была использована технология OpenMP.

Образцы исследования

В качестве образцов исследования используются стандартные для тестирования различных процессов два типа песчаников с сильно различающимися фильтрационно-емкостными свойствами: песчаники Береа и LV60A. Цифровая модель данных пористых сред была заимствована из открытого атласа Imperial College London [12, 13].

Разрешение микротомографического изображения песчаника Береа (рис. 1) составляет 3,2 мкм, а LV60A (рис. 2) — 10 мкм. Числовые значения на координатных осях указаны в микрометрах. Для проведения вычислительных экспериментов из них были вырезаны фрагменты размерами 200x200x200 вокселей.

Рис. 3 — Сравнение кривых капиллярного давления дренирования песчаника Береа, измеренных экспериментально и вычисленных на кубическом фрагменте цифрового керна Fig. 3 — Comparison of drainage capillary pressure curves of Berea sandstone, obtained through experiment and computed on a cube-shaped fragment of digital core

¿¿W -2000 -

18001600-

■ J7M мери МСН i * Молсль

140012001000-еоо-600 -

__

4002000-

■ ■ ■ ■ ■ —■— —1

0,2

0.3

0,4 0.5 0.6 0,7

Насыщсшюсть смачивающей ф;: шп. отн, ил.

Рис. 4 — Сравнение кривых капиллярного давления дренирования песчаника LV60A, измеренных экспериментально и вычисленных на кубическом фрагменте цифрового керна Fig. 4 — Comparison of drainage capillary pressure curves of LV60A sandstone, obtained through experiment and computed on a cube-shaped fragment of digital core

Рис. 5 — Визуализация несмачивающей фазы при дренировании в песчанике Береа. Белой стрелкой отмечены капиллярные пальцы, распространяющиеся в направлениях,

отличных от гидродинамического перепада давления Fig. 5 — Visualization of non-wetting phase while draining in Berea sandstone. White arrows show capillary propagation in directions different from flowing pressure gradient

Пористость песчаников Береа и [УбОА составляет 19,6% и 37% соответственно. Тензор абсолютной проницаемости образцов является практически изотропным, т.к. доля его изотропной компоненты для песчаника Береа 0,670 мкм2, а для [_\/60А — 33 мкм2. Для песчаника Береа компоненты тензора в главных осях следующие: кХХ = 0,69 мкм2, к„ = 0,66 мкм2, кп = 0,66 мкм2; для образца [_\/60А: кж = 34,65 мкм2, куу = 32,75 мкм2, к2 = 32,50 мкм2. На основании методики, описанной в работах [1, 2], было выявлено, что данные объемы образцов являются представительными с точки зрения измерения их фильтрационно-емкостных свойств (ФЕС). Это позволит провести корректное сопоставление вычисленных и экспериментальных кривых капиллярного давления дренирования.

Результаты исследований

Для верификации предложенного в работе метода было выполнено сравнение кривых капиллярного давления, вычисление которых производилось при дренировании на кубических фрагментах керна, с результатами лабораторных исследований методом ртутной порометрии (по поверхностному натяжению). Результаты экспериментов для песчаников Береа и [_\60А представлены в работах [4, 5]. ФЕС цифровых образов песчаников, исследуемых в данной работе, идентичны характеристикам образцов, подвергнутых лабораторным испытаниям. Поверхностное натяжение в расчетах и в эксперименте для песчаника Береа составляло 28,8 мН/м [4], а для [У60А — 50,8 мН/м [5]. Для каждой кривой было вычислено среднее от модуля относительного отклонения по следующей

формуле:

100 Д = N'i

эксперимент

62

где N — число замеров давления (в нашем случае их было 11); 5,:™д™>и 5?1""*""""' — значения насыщенности образца смачивающей фазой при I -ом капиллярном давлении, измеренные в вычислительном и лабораторном экспериментах соответственно.

На рис. 3 и рис. 4 показаны кривые капиллярного давления, полученные экспериментальным методом, и вычисленные по цифровым моделям для двух образцов. Согласно приведенным на данных рисунках результатам, для песчаников Береа и [У60А Л = 8,5% и Л = 3,83% соответственно. Такое относительное отклонение можно считать вполне удовлетворительным.

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

Работа выполнена при поддержке гранта РФФИ мол_а №16-35-00155.

Экспозиция НЕфть газ ноябрь 7 (60) 2017

Итоги

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

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

Выводы

В работе представлена новая методика измерения кривых капиллярного давления по цифровым микротомографическим изображениям. Проведено успешное сопоставление кривых дренирования капиллярного давления, рассчитанных при проведении вычислительного эксперимента на цифровых моделях песчаников Береа и 1_\/60А, и полученных в лабораторных фильтрационных экспериментах.

ENGLISH

OIL PRODUCTION

Measuring of drainage capillary pressure curves with digital microtomography imaging of porous

UDC 622.276

Authors:

Timur R. Zakirov — Ph.D., scientist; timurz0@mail.ru

Akhmet A. Galeev — Ph.D., assistant professor; akhmet.galeev@kpfu.ru

Kazan Federal University, Kazan, Russian Federation

Abstract

The paper explores the new technique of measuring drainage capillary pressure curves using digital microtomography imaging of porous media. Analysis of drainage curves is based on simulation experiments of displacement of one fluid in the pore volume with the other. Computational grid is represented by binary microtomography imaging of sandstone. Mathematical tool used for computation is the lattice Boltzmann method. Calculated curves have been successfully correlated with experimental data. It has been shown that the flow patterns, used in simulation experiments, have revealed peculiar propagation in directions different from flowing pressure gradient, including opposite directions.

Materials and methods

Math modelling of two-phase flows in pore channels incorporates the lattice Boltzmann method. The two-phase boundary line phenomena are described with numeric field gradient. Test samples are represented by digital models of typical Berea sandstone and LV60A sandpack.

Results

The paper presents the new technique of measuring capillary pressure curves using digital microtomography imaging. Drainage capillary pressure curves have been successfully correlated. The curves have been derived from simulation experiments through digital models of Berea sandstone and LV60A sandpack as

well as from filtration experiments in a laboratory.

Conclusions

The paper presents the new technique of measuring capillary pressure curves using digital microtomography imaging. Drainage capillary pressure curves have been successfully correlated. The curves have been derived from simulation experiments through digital models of Berea sandstone and LV60A sandpack as well as from filtration experiments in a laboratory.

Keywords

math modelling, capillary pressure curves, x-ray microtomography, porous media, lattice Boltzmann method

References

1. Mostaghimi P, Blunt M. J., Bijeljic B. Computations of absolute permeability on micro-CT images. Mathematical Geosciences, 2013, Vol. 45, issue 1, pp. 103-125.

2. Zakirov T. R., Galeev A. A., Korolev E. A., Statsenko E. O. Flow properties of sandstone and carbonate rocks by X-ray computed tomography. Current science, 2016, Vol. 110, issue 11, pp. 2142-2147.

3. Blunt M. J., Bijeljic B., Dong H., Gharbi O., Iglauer S., Mostaghimi P., Paluszny A., Pentland C. Pore-scale imaging and modeling. International Journal of Greenhouse Gas Control, 2014, Vol. 22, pp. 1-14.

4. Tsuji T., Jiang F., Christensen K. T. Characterization of immiscible fluid displacement processes with various capillary numbers and viscosity ratios in 3D natural sandstone. Advances in Water Recources, 2016, Vol. 95, pp. 3-15.

5. Ruspini L. C., Farokhpoor R., 0ren P.E. Pore-scale modeling of capillary trapping in water-wet porous media: A new cooperative pore-body filling model. Advances in Water Recources, 2017, Vol. 108, pp. 1-14.

6. Chen S., Doolen G. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics, 1998, V. 30, pp. 329-364.

7. Succi S. Lattice Boltzmann Equation for Fluid Dynamics and Beyond. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford University Press, 2001, 308 p.

8. Reis T., Phillips T.N. Lattice Boltzmann model for simulating immiscible two-phase flows. Journal of Physics A: Mathematical and Theoretical, 2007, V. 40, pp. 4033-4053.

9. Aslan E., Taymaz I., Benim A.C. Investigation of the Lattice Boltzmann SRT and MRT Stability for Lid Driven Cavity Flow. International Journal of Materials, Mechanics and Manufacturing, 2014, Vol. 2, issue 4, pp. 317-324.

10. Huang H., Huang J.-J., Lu X.-Y. Study of immiscible displacements in porous media using a color-gradient-based multiphase lattice Boltzmann method. Computers & Fluids, 2014, issue 93, pp, 164-172.

11. Zou Q, He X. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Physics of Fluids, 1997, Vol. 9, issue 6, pp. 1591-1598.

12. Imperial College London. Berea sandstone. Available at: http://www. imperial.ac.uk/earth-science/research/ research-groups/perm/research/ pore-scale-modelling/micro-ct-images-and-networks/berea-sandstone/

13. Imperial College London. Sand pack LV60A. Available at: http:// www.imperial.ac.uk/ earth-science/research/research-groups/perm/research/pore-scale-modelling/micro-ct-images-and-networks/sand-pack-lv60a/