CC BY
12УДК 531.724:004.941+622.276:519.8
МОДЕЛИРОВАНИЕ И ВИЗУАЛИЗАЦИЯ ПРОЦЕССА ВЫТЕСНЕНИЯ НЕФТИ ИЗ ПОРИСТОЙ СРЕДЫ
М. В. Михайлюк, П. Ю. Тимохин, А. В. Мальцев, В. Ф. Никитин, Е. И. Скрылева,
В. В. Тюренкова*
Федеральный научный центр Научно-исследовательский институт системных исследований
Российской академии наук *tyurenkova. v. v@yandex. ru
В данной работе предлагаются методы моделирования и визуализации процесса вытеснения вязкой жидкости из образца пористой среды. Моделирование выполняется на основе специальной системы уравнений, описывающих взаимодействие в пористой среде двух жидкостей, имеющих различную вязкость. Рендеринг полигональной модели изоповерхности насыщенности вытесняющей жидкости осуществляется в масштабе реального времени с использованием возможностей современных графических карт и разработанной собственной модификации метода «шагающих кубиков» (marching cubes). На основе предлагаемых решений реализован программный комплекс с эргономичным пользовательским интерфейсом.
Ключевые слова: вытеснение, фильтрация, визуализация, поверхность уровня, OpenGL.
MODELING AND VISUALIZATION OF PROCESS OF OIL DISPLACEMENT FROM POROUS MEDIUM
M. V. Mikhaylyuk, P. Yu. Timokhin, A. V. Maltsev, V. F. Nikitin, E. I. Skryleva,
V. V. Tyurenkova*
Institute for System Analysis, Russian Academy of Sciences, * tyurenkova.v.v@yandex.ru
The paper suggests simulation and visualization methods of the viscous fluid displacement process from a porous medium sample. The simulation is performed on the basis of special system of equations describing the interaction of two liquids with different viscosities in a porous medium. A polygonal model of the displacing fluid's saturation iso-surface is rendered in real-time using the capabilities of modern graphics cards and an original modification of the marching cubes method. A user-friendly software product implementing the proposed solutions has been developed.
Keywords: displacement, filtration, visualization, level surface, OpenGL.
Введение
Фильтрационные течения вязких жидкостей в пористых средах являются определяющими в процессах добычи жидких полезных ископаемых (нефти, газоконденсата и пр.). Повышение нефтеотдачи нефтеносных пластов связано с качеством вытеснения нефти, насыщающей подземные резервуары, с помощью воды или полимерных растворов.
Известно, что вытеснение более вязкой жидкости менее вязкой жидкостью из пористой среды является неустойчивым процессом. Возникающая неустойчивость приводит к нарушению формы поверхности раздела фаз и переходу от фронтального вытеснения к прорыву отдельных языков вытесняющей жидкости и захвату вытесняемой жидкости внутри пористой среды. Эти эффекты существенно снижают качество вытеснения. Присутствующие изначально небольшие неоднородности проницаемости в пористой среде при вынужденных интенсивных фильтрационных течениях могут также приводить к формированию трещин и каналов в пористой матрице, что увеличивает неоднородность проницаемости [1, 2].
При моделировании описанного процесса одной из важных задач является визуализация полученных результатов, в частности, поверхности раздела фаз. Возникающая неустойчивость приводит к образованию «пальцев» различной длины и толщины, которые
могут как разветвляться, так и сливаться в один. Также возможно образование объемов жидкости, которые не поддаются вытеснению. Реализация трехмерной модели поверхности раздела фаз имеет ряд сложностей, преодоление которых позволит более точно и качественно описать процесс фильтрации в пористой среде.
3D визуализация дает возможность построить полигональную модель поверхности уровня (изоповерхности) изучаемой величины и синтезировать ее изображение в масштабе реального времени. При этом обеспечивается возможность on-line перемещения пользователем виртуальной камеры с целью изучения поверхности с разных сторон.
Для моделирования и визуализации поверхности уровня можно использовать метод «бросания лучей» (ray casting) [3], который позволяет изучать внешнюю (видимую) поверхность уровня. В настоящей работе предлагается использовать метод «шагающих кубиков» (marching cubes), строящий также и внутренние поверхности уровня, которые можно анализировать с использованием сечений (например, плоскостью).
Входное сечение
-Н/2
Моделирование процесса вытеснения
Предполагается, что область пористой среды имеет форму прямоугольного параллелепипеда, заполненного вытесняемой жидкостью (рис.1). Через одно из сечений подается жидкость меньшей вязкости. Выход осуществляется через противоположную грань области. На остальных гранях задано условие непротекания.
Пористость, проницаемость, вязкость каждой жидкости, поверхностное натяжение и угол смачивания считаются постоянными. Также предполагается, что скорость жидкости достаточно мала, температура постоянна и межфазного массообмена нет. Капиллярное давление считается по модели Брукса-Кори р^ - СБ—, где константы С и а определяются эмпирическим путем [4].
Определяющая система уравнений для данной задачи, описанная подробно в статье [5], выгладит следующим образом.
Закон сохранения массы для каждой жидкости: р, д ,
-н--(ф^рик ) - 0, где / -1,2 соответствует вытесняемой и вытесняющей
дг дх. ,;
жидкости, ф - пористость, ^ - насыщенность, р - плотность (истинная), и. - компонента у истинной скорости фазы;
Закон Дарси для каждой фазы, учитывающий капиллярной давление рС :
Рис. 1
Выходное сечение
<$sk Pk =
К KR dpk
^ дХу
Здесь ^ - вязкость фазы, К£ - относительная проницаемость, К0 - абсолютная проницаемость среды, р^ - внутрипоровое давление в фазе.
Разность фаз равна капиллярному давлению р - р2 - рС (.
Полученные уравнения решаются методом конечных разностей с размещением всех параметров в узлах трехмерной равномерной сетки, каждый распределенный параметр представляет собой трехмерный массив. Расчет на каждом шаге по времени проводится явно-неявными методами, в том числе Bi-CGStab для уравнения, рассчитывающего давление фаз (неявного) и метода TVD для продвижения насыщенности фаз во времени (явного).
После каждого шага формируются массивы давления и насыщенности вытесняющей фазы, которые далее подвергаются обработке, как для визуализации результатов, так и для расчета интегральных параметров процесса.
Метод визуализации результатов моделирования
Исходными данными для визуализации являются последовательность трехмерных массивов вещественных чисел, полученных в результате расчета математической модели, описанной в предыдущем разделе. Каждый массив записан в отдельном бинарном файле и соответствует определенному временному шагу расчета математической модели. Размеры массива совпадают с размерами расчетной сетки, а значениями элементов массива являются насыщенности ^ к вытесняющей жидкости в ячейках (i, j, к) расчетной сетки.
Визуализация результатов моделирования выполняется в разработанной программной оболочке. Оболочка с определенным временным интервалом загружает в оперативную память каждый следующий трехмерный массив насыщенностей, рассчитывает на его основе полигональную модель изоповерхности для заданного пользователем уровня s = s* насыщенности и визуализирует полученную модель в отдельном окне.
В окне визуализации пользователь имеет возможность интерактивно взаимодействовать с полигональной моделью изоповерхности (произвольно вращать, приближать/удалять и др.). Во избежание возникновения задержек при взаимодействии пользователя с моделью изоповерхности, загрузка массивов насыщенностей и расчет полигональных моделей изоповерхностей выполняются в отдельном потоке.
Для расчета полигональной модели изоповерхности в данной работе предложена реализация алгоритма «шагающих кубиков», обеспечивающая совпадение общих вершин смежных треугольников в результирующей модели.
Построение полигональной модели изоповерхности
Опишем кратко метод marching cubes построения изоповерхности трехмерного скалярного поля (более полное описание см. в [6, 7]). В объеме моделирования определим новую равномерную сетку с размером ячейки R х R х R, где R - некоторый вещественный параметр. Будем называть эту сетку сеткой рендеринга. Ячейку данной сетки будем идентифицировать индексами (i, j, к) , вершины каждой ячейки пронумеруем целыми числами от 0 до 7, а ребра - от 0 до 11. Тем вершинам сетки рендеринга, для которых значения соответствующих им ячеек моделирования превышают s*, припишем битовое значение bF равное 1. Оставшимся вершинам - значение Ък = 0. Тогда любую конфигурацию из 0 и 1, приписанных вершинам каждой ячейки, можно представить в виде 8 битов, т. е. одного байта, который мы будем обозначать KtJ к. Каждой такой конфигурации
ячейки сетки рендеринга однозначно соответствует набор треугольников (полигонов), который будет являться частью строящейся изоповерхности. Этот набор запишем в массив T как элемент с номером K . к . Вершины полигонов набора находятся только на тех ребрах
ячеек сетки, на концах V и V2 которых bv Ф bv, при этом положение этих вершин
определяется с помощью линейной интерполяции на основе значений s*, s(V ) и s(V2 ). Алгоритм построения полигонов таков, что на каждом ребре ячейки будет не более одной вершины полигона, и полигоны соседних ячеек будут иметь общие вершины или ребра. Тем самым, полученная в результате полигональная модель будет замкнутой и не содержащей дыр. При этом внутри одной поверхности уровня могут находиться одна или несколько других поверхностей уровня. Их можно будет увидеть только при реализации сечений модели (например, плоскостью). Для реализации алгоритма необходимы следующие структуры данных:
• Массив МЕ, элементы которого соответствуют ребрам сетки рендеринга.
Первоначально он пустой, а затем в те элементы, которые соответствуют ребрам, содержащим вершину полигона, запишутся координаты этой вершины и нормаль N в ней.
• Массив Т из 256 элементов, соответствующих всевозможным конфигурациям ячеек. В каждом элементе записаны локальные номера ребер треугольников (от 0 до 11) участка изоповерхности для такой конфигурации.
• Массив С размером равный числу ячеек сетки рендеринга, элементами которого являются значения конфигураций К этих ячеек.
• Массив П всех полигонов изоповерхности, в который будут записываться полигоны участков изоповерхности внутри ячеек сетки рендеринга. Каждый полигон представлен индексами трех вершин - элементов массива М£, причем вершины перечислены в таком порядке, который однозначно определяет направление нормали к полигону, соответствующее его лицевой стороне.
Заметим, что координаты конкретных вершин или ребер сетки рендеринга хранить не надо, т.к. их легко вычислить, зная размер ячейки. Таким образом, мы получаем следующий
Алгоритм построения массива П полигонов изоповерхности
Начало.
1. Создаем пустой массив М£ по числу всех ребер сетки (в порядке перебора индексов (г, у, к)).
2. Для каждой вершины V сетки вычисляем значение ЪГ .
3. Для каждой ячейки С^к сетки рендеринга вычисляем конфигурацию К и записываем в массив С .
4. Цикл по всем ребрам Е . к сетки рендеринга
Если для вершин V и V этого ребра Ъу Ф Ъ^, то вычисляем на ребре точку
- VI н £ - ^1^2 - V).
1 з(У2) -
Записываем координаты р ^ в элемент массива Мя, соответствующий Е^к,
и ставим флаг, что это ребро содержит точку. Конец цикла. 5. Цикл по всем элементам С массива С
Получаем значение К элемента С .
В массиве Т выбираем элемент с индексом К , в котором записаны номера
ребер треугольников (полигонов) участка изоповерхности, попадающего в эту ячейку.
По этим номерам вычисляем соответствующие элементы массива М . Тройки номеров элементов массива М записываем в массив П полигонов изоповерхности.
Для каждого записанного полигона вычисляем нормаль и добавляем ее к нормалям вершин этого полигона в массиве МЕ.
Конец цикла. Конец.
После конца работы алгоритма в заполненных элементах массива МЕ поле N нормали будет содержать сумму нормалей ко всем полигонам, содержащим эту вершину. Естественно, эта сумма в общем случае не будет нормализованной.
Визуализация полигональной модели
Для визуализации построенной полигональной модели поверхности уровня предлагается использовать графическую библиотеку OpenGL. С помощью функций этой библиотеки можно установить и направить в нужную сторону источник освещения, задать свойства этого источника, а также свойства материала, с помощью которого будет моделироваться изоповерхность (см. [8, 9]). Кроме того, можно задать некоторые параметры расчета освещенности.
Установить и направить источник освещения можно с помощью операторов glLightfv(GL_LIGHT0, GL_POSITION, pos), glLightf(GL_LIGHT0, GL_SPOT_DIRECTION, dir), где pos и dir - трехмерные векторы положения и направления соответственно.
Интенсивности испускаемого фонового, диффузного и зеркального излучения, а также сфокусированность источника задаются операторами glLightfv(GL_LIGHTi, GL_AMBIENT, Ia), glLightfv(GL_LIGHTi, GL_DIFFUSE, Id), glLightfv(GL_LIGHTi, GL_SPECULAR, Is), glLightf(GL_LIGHTi, GL_SPOT_EXPONENT, t), Id и Is задают 3D векторы цвета в формате (r, g, b).
Аналогично, с помощью операторов glMaterialfv(GL_FRONT, GL_AMBIENT, ka), glMaterialfv(GL_FRONT, GL_DIFFUSE, kd), C^ç^^ / glMaterialfv(GL_FRONT, GL_SPECULAR, ks),
^S /а У glMaterialfv(GL_FRONT, GL_SPECULAR, spec),
/ yS можно задать коэффициенты ка, kd и ks отражения
фоновой, диффузной и зеркальной составляющей света, попадающего на лицевую поверхность поверхности уровня и степень p зеркального отражения материала.
Тогда освещенность IV в вершине V полигональной сетки вычисляется по формуле
где I
P
ф
Iv = (Iaka + Idkd max(0, Cos6) + Isks max(0, Cosa)p )Cos'ф.
Рис. 2
где 6 - угол между нормалью к поверхности в точке V и направлением u на источник освещения (см. рис. 2), а - угол между отраженным от поверхности лучом r и вектором s, направленным на наблюдателя, ф - угол между направлением m источника освещения и направлением на вершину V.
Для расчета косинуса всех углов достаточно знать координаты расположения наблюдателя и источника освещения. Действительно, тогда можно определить все необходимые векторы п, т, и, s и г, нормализовать их и вычислить скалярные произведения
CosQ = {n,u), Cosa = (r,s), Cosy = (m,u).
Нормализованный вектор Я нормали в точке Р получается из нормали Nр в той же
точке по формуле п = N /
NT
На рис. 3 показан пример визуализации модели изоповерхности насыщенности вытесняющей жидкости для = 0.15 .
Построенную поверхность можно изучать в масштабе реального времени, перемещая виртуальную камеру с помощью клавиатуры. Реализовано перемещение вдоль и повороты камеры вокруг ее осей координат.
Для ускорения работы алгоритма его целесообразно выполнять на графической карте, используя технологию шейдеров.
Заключение
В данной работе предложена эффективная технология моделирования и визуализации процесса вытеснения вязкой жидкости из образца пористой среды, позволяющая синтезировать в масштабе реального времени динамическую интерактивную 3D-модель изоповерхности насыщенности вытесняющей жидкости с применением современных технологий обработки данных на графических ускорителях. В этой технологии моделирование процесса вытеснения выполняется на основе специальной системы уравнений, описывающих взаимодействие в пористой среде двух жидкостей, имеющих различную вязкость. На каждом шаге моделирования в результате рассчитывается 3D-массив насыщенностей вытесняющей жидкости, на основе которого строится полигональная модель изоповерхности. Это осуществляется с помощью предложенной в данной работе модификации метода «шагающих кубиков», позволяющей строить замкнутые полигональные модели изоповерхностей, не содержащие дыр. Были созданы методы и алгоритмы хранения, обработки и визуализации синтезируемой полигональной модели Рис. 3
изоповерхности, работающие в масштабе реального
времени и использующие графическую библиотеку OpenGL. Также были предложены методы и алгоритмы, основанные на распределенной обработке потоков данных моделирования и визуализации, позволяющие интерактивно взаимодействовать с синтезированной 3D-моделью изоповерхности (произвольно вращать/приближать, изменять значение изоуровня и др.) и анимировать динамику ее изменения в ходе моделирования процесса вытеснения. Созданная технология, методы и алгоритмы реализованы в виде программного комплекса и успешно апробированы на ряде сложных конфигураций реализованной математической модели вытеснения. Предложенное решение позволяет заметно повысить эффективность анализа результатов математического моделирования процесса вытеснения, а также оперативно разрабатывать и исследовать новые конфигурации реализованной математической модели вытеснения в направлении повышения нефтеотдачи нефтеносных пластов.
Работа выполнена при поддержке гранта РФФИ № 16-29-15099 офи_м.
Литература
1. Smirnov N. N., Legros J. C., Nikitin V. F., Istasse E., Schramm L., Wassmuth F., Hart D. A. Filtration in artificial porous media and natural sands under microgravity conditions// Microgravity Science and Technology, (2003) 14 (2), pp. 3-28.
2. Smirnov N. N., Nikitin V. F., Ivashnyov O. E., Maximenko A., Thiercelin M., Vedernikov A., Scheid B., Legros J. C. Microgravity investigations of instability and mixing flux in frontal displacement of fluids // Microgravity Science and Technology, (2004) 15 (2), pp. 35-51.
3. Смирнов Н. Н., Никитин В. Ф., Михайлюк М. В., Тимохин П. Ю., Тюренкова В. В., Стамов Л.И. Визуализация результатов моделирования неустойчивого вытеснения нефти из пористых сред // Труды НИИСИ РАН. 2016. Т.6. №2.
4. Smirnov N.N., Nikitin V.F., Norkin A.V., Kiselev A.B., Legros J.C., Istasse E. "Microgravity investigation of capillary forces in porous media" // Space Forum. 2000. № 6 (1-4). pp. 1-10.
5. Козлов И. В., Скрылева Е. И. Математическое моделирование и обработка эксперимента по вытеснению нефти водой из неокомских песчаников // Вестник кибернетики. 2016. №2. С. 138-145.
6. Marching cubes. URL: https://ru.wikipedia.org/wiki/Marching_cubes (дата обращения: 02.08.2016).
7. Polygonising a scalar field. URL: http://paulbourke.net/geometry/polygonise/ (дата обращения: 02.08.2016).
8. Михайлюк М. В. Основы компьютерной графики. М.: МГТУ МИРЭА. 2011. - 156 с.
9. Хилл Ф. OpenGL. Программирование компьютерной графики. - СПб.: Питер, 2002. - 1088 с.
