МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ФИЛЬТРАЦИИ В ПЕРКОЛЯЦИОННЫХ РЕШЕТКАХ С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ СВЕРХВЫСОКОЙ ПРОИЗВОДИТЕЛЬНОСТИ
52 М. Н. Воронюк
УДК 519.6
М. Н. ВОРОНЮК
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ФИЛЬТРАЦИИ В ПЕРКОЛЯЦИОННЫХ РЕШЕТКАХ
С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ СВЕРХВЫСОКОЙ ПРОИЗВОДИТЕЛЬНОСТИ
Рассматривается процесс математического моделирования заводнения нефтяного месторождения в рамках динамической перколяционной модели на решетках, содержащих 109 узлов, с использованием суперкомпьютеров с GPU. Предложен алгоритм эффективной балансировки загрузки вычислительных узлов.
Ключевые слова: нефтедобыча, GPU, перколяция, псевдослучайные числа, метод Монте-Карло, балансировка загрузки.
Введение. Современные методы добычи углеводородов невозможны без использования
высокопроизводительных автоматизированных систем управления. Такие автоматизирован-
ные системы, базирующиеся на статистических моделях, позволяют оценить связность двух и
более скважин и проводимость между ними, время прорыва агента между добывающей и на-
гнетательной скважинами и его поведение после прорыва (например, время достижения
50%-ной обводненности вокруг добывающей скважины) [1].
Использование данных моделей сопряжено, однако, с большими вычислительными за-
тратами. Применение для моделирования вычислительных систем сверхвысокой производи-
тельности, содержащих графические ускорители (Graphics Processing Unit — GPU), позволяет
за приемлемое время получать результаты с высокой точностью, но для эффективного прове-
дения расчетов на графических ускорителях используемые алгоритмы требуют адаптации.
Постановка задачи. Для описания процессов нефтевытеснения используется динами-
ческая перколяционная модель (ДПМ) [2]. Пласт представляется в виде агрегата из трех взаи-
мопроникающих составляющих: скелета, жидкой фазы и газовой фазы. В модели нефтесо-
держащие поры и капиллярные каналы представлены в виде правильной кубической решетки, содержащей 109 узлов. Каждый канал случайным образом помечается „открытым“ с веро-
ятностью протекания жидкости p (0≤p≤1) или „закрытым“ с вероятностью (1–p). Открытым
каналам приписывается значение гидросопротивления. Состояние (α) каждого узла (поры)
описывается линейной комбинацией трех состояний: узел заполнен нефтью: α= +1, узел за-
полнен водой: α = –1 или узел не заполнен (вакантен): α = 0.
Процесс фильтрации описывается системой кинетических уравнений для решеточных
функций [2]:
∑d
dt
Gm
(t)
=
−
γm
ηγm Ωγm
(Gm
−
Gm + γ m
)
;
∑d
dt
Fmα (t)
=
−
γm
ηγm Ωγm
(Ψm+γm ,m Fmα+γm Fm0
−
Ψ m ,m + γ m
Fm0+γm Fmα )
;
∑Fm0 (t) = 1 −
Fmα (t) ,
α=±1
где m — векторный индекс с целочисленными координатами, характеризующий узел с ин-
дексом ijk; величины Fmα (t) , Fm0 (t) интерпретируются как средние доли нефти, воды и вакан-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 5
Математическое моделирование процессов фильтрации
53
сии в узле m; Gm (t) — давление в узле m; γm — вектор, нумерующий ребра, выходящие из
узла m (одна из координат вектора γm равна ±1, остальные — нулевые);
( )Ψm+γm,m = Gm+γm − Gm − (Gm+γm − Gm ) 2 — функция-„ниппель“ ребра γm; ηγm — случай-
ная величина со значениями „0“ и „1“ представляет собой индикатор открытости канала γm;
Ωγm — безразмерная случайная проницаемость канала, нормированная к единице.
Для численного моделирования процессов фильтрации использована следующая аппроксимация:
∑Gm = Gm − τ ηγm Ωγm (Gm − Gm+γm ) ; γm
(1)
∑Fmα = Fmα − τ ηγm Ωγm (Ψm+γm ,m Fmα+γm Fm0 − Ψm,m+γm Fm0+γm Fmα ) ; γm
(2)
∑Fm0 = 1−
Fmα (t) ,
α=±1
(3)
где τ — шаг по времени; Gm , Fmα , Fm0 и Gm , Fmα , Fm0 — значения решеточной функции в
моменты времени t и t+τ соответственно.
Методика расчета. Рассматриваемая задача (1)—(3) является стохастической и требует
проведения многократных расчетов каждого из сценариев заводнения нефтяного месторож-
дения. Соответственно необходимо задание нового набора случайных величин, характери-
зующих гидросопротивление межпоровых каналов, для каждого из сценариев, и, кроме того,
само их моделирование требует значительных затрат процессорного времени.
Использование внутреннего параллелизма, присущего динамической перколяционной
модели, позволяет эффективно производить расчеты на основе параллельных систем сверх-
высокой производительности.
Перед непосредственным расчетом сценариев заводнения выполняются следующие под-
готовительные действия, что позволяет, с одной стороны, избежать их многократного повто-
рения, а с другой — эффективно использовать параллельные системы:
1) генерация геометрии месторождения с помощью механизма R-функций [3];
2) „наложение“ скважин;
3) решение задачи балансировки загрузки вычислительных узлов.
Информация о геометрии месторождения задается с помощью базовых элементов раз-
личных форм (эллипсоида, параллелепипеда, эллиптического цилиндра, полупространства) и
операций над ними (пересечение, объединение, дополнение, поворот, сдвиг).
Скважины задаются набором цилиндрических секций заданной толщины.
При окончательном расчете каждого из сценариев заводнения выполняются следующие
операции:
1) генерация перколяционной решетки с помощью генератора псевдослучайных чисел;
2) „наложение“ скважин на перколяционную решетку;
3) решение системы дифференциальных уравнений (1)—(3).
При генерации перколяционной решетки для задания значения случайной величины
ηγm используется генератор псевдослучайных чисел на основе М-последовательностей [4, 5]
ϕ j+1023 = ϕ j+511 + ϕ j+127 + ϕ j+7 + ϕ j mod 2 ,
(4)
где ϕi — одноразрядные члены рекуррентной последовательности, и применяется формула для получения очередного псевдослучайного 32-битового числа [5]:
( ) ( ) ( ) ( )wj+8 = wj > 31 ⊕ wj > 24 ⊕ wj+4 ⊕ wj+16 ,
(5)
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 5
54 М. Н. Воронюк
где ⊕ — операция побитового сложения по модулю 2 целых чисел; операции „>>“ и „
УДК 519.6
М. Н. ВОРОНЮК
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ФИЛЬТРАЦИИ В ПЕРКОЛЯЦИОННЫХ РЕШЕТКАХ
С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ СВЕРХВЫСОКОЙ ПРОИЗВОДИТЕЛЬНОСТИ
Рассматривается процесс математического моделирования заводнения нефтяного месторождения в рамках динамической перколяционной модели на решетках, содержащих 109 узлов, с использованием суперкомпьютеров с GPU. Предложен алгоритм эффективной балансировки загрузки вычислительных узлов.
Ключевые слова: нефтедобыча, GPU, перколяция, псевдослучайные числа, метод Монте-Карло, балансировка загрузки.
Введение. Современные методы добычи углеводородов невозможны без использования
высокопроизводительных автоматизированных систем управления. Такие автоматизирован-
ные системы, базирующиеся на статистических моделях, позволяют оценить связность двух и
более скважин и проводимость между ними, время прорыва агента между добывающей и на-
гнетательной скважинами и его поведение после прорыва (например, время достижения
50%-ной обводненности вокруг добывающей скважины) [1].
Использование данных моделей сопряжено, однако, с большими вычислительными за-
тратами. Применение для моделирования вычислительных систем сверхвысокой производи-
тельности, содержащих графические ускорители (Graphics Processing Unit — GPU), позволяет
за приемлемое время получать результаты с высокой точностью, но для эффективного прове-
дения расчетов на графических ускорителях используемые алгоритмы требуют адаптации.
Постановка задачи. Для описания процессов нефтевытеснения используется динами-
ческая перколяционная модель (ДПМ) [2]. Пласт представляется в виде агрегата из трех взаи-
мопроникающих составляющих: скелета, жидкой фазы и газовой фазы. В модели нефтесо-
держащие поры и капиллярные каналы представлены в виде правильной кубической решетки, содержащей 109 узлов. Каждый канал случайным образом помечается „открытым“ с веро-
ятностью протекания жидкости p (0≤p≤1) или „закрытым“ с вероятностью (1–p). Открытым
каналам приписывается значение гидросопротивления. Состояние (α) каждого узла (поры)
описывается линейной комбинацией трех состояний: узел заполнен нефтью: α= +1, узел за-
полнен водой: α = –1 или узел не заполнен (вакантен): α = 0.
Процесс фильтрации описывается системой кинетических уравнений для решеточных
функций [2]:
∑d
dt
Gm
(t)
=
−
γm
ηγm Ωγm
(Gm
−
Gm + γ m
)
;
∑d
dt
Fmα (t)
=
−
γm
ηγm Ωγm
(Ψm+γm ,m Fmα+γm Fm0
−
Ψ m ,m + γ m
Fm0+γm Fmα )
;
∑Fm0 (t) = 1 −
Fmα (t) ,
α=±1
где m — векторный индекс с целочисленными координатами, характеризующий узел с ин-
дексом ijk; величины Fmα (t) , Fm0 (t) интерпретируются как средние доли нефти, воды и вакан-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 5
Математическое моделирование процессов фильтрации
53
сии в узле m; Gm (t) — давление в узле m; γm — вектор, нумерующий ребра, выходящие из
узла m (одна из координат вектора γm равна ±1, остальные — нулевые);
( )Ψm+γm,m = Gm+γm − Gm − (Gm+γm − Gm ) 2 — функция-„ниппель“ ребра γm; ηγm — случай-
ная величина со значениями „0“ и „1“ представляет собой индикатор открытости канала γm;
Ωγm — безразмерная случайная проницаемость канала, нормированная к единице.
Для численного моделирования процессов фильтрации использована следующая аппроксимация:
∑Gm = Gm − τ ηγm Ωγm (Gm − Gm+γm ) ; γm
(1)
∑Fmα = Fmα − τ ηγm Ωγm (Ψm+γm ,m Fmα+γm Fm0 − Ψm,m+γm Fm0+γm Fmα ) ; γm
(2)
∑Fm0 = 1−
Fmα (t) ,
α=±1
(3)
где τ — шаг по времени; Gm , Fmα , Fm0 и Gm , Fmα , Fm0 — значения решеточной функции в
моменты времени t и t+τ соответственно.
Методика расчета. Рассматриваемая задача (1)—(3) является стохастической и требует
проведения многократных расчетов каждого из сценариев заводнения нефтяного месторож-
дения. Соответственно необходимо задание нового набора случайных величин, характери-
зующих гидросопротивление межпоровых каналов, для каждого из сценариев, и, кроме того,
само их моделирование требует значительных затрат процессорного времени.
Использование внутреннего параллелизма, присущего динамической перколяционной
модели, позволяет эффективно производить расчеты на основе параллельных систем сверх-
высокой производительности.
Перед непосредственным расчетом сценариев заводнения выполняются следующие под-
готовительные действия, что позволяет, с одной стороны, избежать их многократного повто-
рения, а с другой — эффективно использовать параллельные системы:
1) генерация геометрии месторождения с помощью механизма R-функций [3];
2) „наложение“ скважин;
3) решение задачи балансировки загрузки вычислительных узлов.
Информация о геометрии месторождения задается с помощью базовых элементов раз-
личных форм (эллипсоида, параллелепипеда, эллиптического цилиндра, полупространства) и
операций над ними (пересечение, объединение, дополнение, поворот, сдвиг).
Скважины задаются набором цилиндрических секций заданной толщины.
При окончательном расчете каждого из сценариев заводнения выполняются следующие
операции:
1) генерация перколяционной решетки с помощью генератора псевдослучайных чисел;
2) „наложение“ скважин на перколяционную решетку;
3) решение системы дифференциальных уравнений (1)—(3).
При генерации перколяционной решетки для задания значения случайной величины
ηγm используется генератор псевдослучайных чисел на основе М-последовательностей [4, 5]
ϕ j+1023 = ϕ j+511 + ϕ j+127 + ϕ j+7 + ϕ j mod 2 ,
(4)
где ϕi — одноразрядные члены рекуррентной последовательности, и применяется формула для получения очередного псевдослучайного 32-битового числа [5]:
( ) ( ) ( ) ( )wj+8 = wj > 31 ⊕ wj > 24 ⊕ wj+4 ⊕ wj+16 ,
(5)
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 5
54 М. Н. Воронюк
где ⊕ — операция побитового сложения по модулю 2 целых чисел; операции „>>“ и „