Параллельная численная реализация пространственно- трехмерной модели фильтрации двухфазной жидкости на основе усовершенствованного попеременно-треугольного метода*

Abstract

В рамках модели Баклея-Леверетта рассмотрена постановка пространственнотрехмерной задачи фильтрации двухфазной несжимаемой жидкости с учетом силы тяжести и без учета капиллярных сил в переменных «давление-водонасыщенность». Для дискретной модели, аппроксимирующей гранично-начальную задачу фильтрации, рассмотрен усовершенствованный итерационный модифицированный попеременно-треугольный метод (МПТМ), учитывающий наличие источников (стоков) в относительно малом числе узлов, имеющих характер дельта-образных (обобщенных) функций. Выполнено численное решение вычислительно трудоемкой задачи в пласте с существенной неоднородностью проницаемости, изменяющейся на четыре порядка в пределах пласта, который имеет значительное превосходство по требуемому для достижения заданной точности числу итераций - в 10 и более раз на сетках, включающих до 3200000 узлов. Для разработанного последовательного варианта, усовершенствованного МПТМ, предложен параллельный алгоритм его численной реализации на многопроцессорной системе с распределенной памятью - вычислительной системе Южного федерального университета в г. Таганроге, содержащей 2048 ядер, c пиковой производительностью 18,8 Тфлопс, базирующийся на методе декомпозиции по двум пространственным направлениям. Показано существенное преимущество предложенного параллельного алгоритма, по сравнению с традиционноиспользуемыми для численного решения данных задач - время моделирования процесса заводнения пластовой системы сократилось в десятки раз. Ключевые слова: задача фильтрации, водонасыщенность, модифицированный попеременно-треугольный метод Математическое описание процессов фильтрации при проектировании разработок нефтяных месторождений приводит к необходимости решения начально- краевых задач для нелинейных уравнений с частными производными, которые в случае трехмерной постановки задачи двухфазной фильтрации (вода-нефть) без учета сжимаемости отдельных фаз, могут быть записаны в виде [1]:  f1(s)kh / 1  f2 (s)kh /  2  px x   f1(s)kh / 1  f2 (s)kh /  2  py y    f1(s)kv  pz  g1  / 1  f2 (s)kv  pz  g 2  / 2 z  q1  0 , mst   f2 (s)kh px /  2 x   f2 (s)kh py /  2 y   f2 (s)kh  pz  g 2  /  2 z  q2 , где s  s(x, y) - водонасыщенность; p  p(x, y,t) - давление; f1(s), f2 (s) - относительные фазовые проницаемости для нефти и воды соответственно; m - пористость пласта; 1,  2 - вяз* Работа выполнена при частичной поддержке Задания №2014/174 в рамках базовой части государственного задания Минобрнауки России

-

1. Введение ( 1 ) ( 2 ) кость нефти и воды соответственно; k(x, y) – проницаемость пласта; функции q1 и q2 – функции, моделирующие работу скважин; 1 и  2 – плотности нефти и воды соответственно, g – ускорение свободного падения.

Граничные условия будем рассматривать в связи с тем, является ли граница  области G проницаемой или нет. На непроницаемой границе поток по нормали должен быть равен нулю, что приводит к условию pn  0 . На проницаемой границе, которой также является контур, ограничивающий скважину, будем рассматривать граничные условия 1 и 2 рода. При совместном движении фаз, которые мы будем предполагать равными wn  w1  w2 , где w1, w2  потоки нефти и воды соответственно, удовлетворяющие условиям: wi  kfi (s)grad p / i . Граничные условия 2 рода реализуются, когда задан поток, закачиваемый (отбираемый) на части цилиндрической границы 1 :

wn (x, y, z,t)  W (x, y, z,t), s(x, y, z,t)  S(x, y, z,t), (x, y, z)  1, 0  t  t0. а также для заданного отбора, заданного на части цилиндрической границы 2 : wn (x, y,t)  W (x, y,t), (x, y, z)  2, 0  t  t0, или давления: wn (x, y,t)  W (x, y,t), p(x, y,t)  P(x, y,t), (x, y, z)  2 , 0  t  t0.

Разумеется, на части цилиндрической поверхности 3 также могут быть заданы условия непроницаемости. Будем использовать при построении дискретных аналогов уравнений ( 1 ) и на скважинах задавать либо дебиты, либо забойные давления. В качестве функциональных зависимостей для задания fi (s), i  1, 2 будем использовать полиномы второго-третьего порядков. 2. Явно-неявная дискретная модель

Будем считать, что G(x, y, z)  область, в которой находится решение начально-граничной задачи для уравнений ( 1 ),( 2 ) – параллелепипед. В области G используем равномерную пространственную сетку h :

h   xi  ihx , y j  jhy , zk  khz ,0  i  Nx , 0  j  N y ,0  k  Nz  ( 3 ) и неравномерную временную сетку ˆ . Разностные аппроксимации уравнений ( 1 ) и ( 2 ) во  внутренних узлах сетки h имеют вид:

B( 1 ) i1/2, j,k

Pi1, j,khx2 Pi, j,k  Bi11/2, j,k Pi, j,k hx2Pi1, j,k  Bi(,2j)1/2,k Pi, j1,khy2 Pi, j,k  Bi(,2j)1/2,k Pi, j1,khy2 Pi, j,k  Bi(,3j),k1/2 Pi, j,k1hz2 Pi, j,k  Bi,3j,k1/2 Pi, j,k hz2Pi, j,k1    g1 / 1   f1,i, j,k1  f1,i, j,k  kv ,i, j,k 1/2   f1,i, j,k 1  f1,i, j,k  kv, i, j,k 1/2 

2hz   g 2 /  2   f2,i, j,k1  f2,i, j,k  kv,i, j,k 1/2   f2,i, j,k 1  f2,i, j,k  kv,i, j,k 1/2  f p,i, j,k ,

2hz где f p  0 вне скважин; f p  Bg  Pw  P на нагнетательных скважинах; f p  Bg  Po  P на эксплуатационных скважинах; Pw , Po – забойные давления на нагнетательной и эксплуатационной скважинах соответственно; Bg   kh hz  f1(s) / 1  f2 (s) / 2  ,   2 / ln h / 2 2rc   c .

Явная схема для аппроксимации уравнения, описывающего изменение водонасыщенности, имеет вид: mH sin,j,1k  si, j,k  1 f2,i1, j,k kh,i1, j,k 1 f2,i, j,k kh,i, j,k  Pi1, j,k 2Pi, j,k 

n  n1  2hx  l  (1  l ) , 1, 2,3, 4,5,6 определяются из условия ориентации схемы против потока.

Соотношения, определяющие коэффициенты разностных схем ( 4 ), ( 5 ), получены интегроинтерполяционным методом и для краткости опущены. 3. Усовершенствованный итерационный МПТМ

Перепишем систему разностных уравнений ( 4 ), ( 5 ) в виде с однородными граничными условиями, изменив соответствующим образом правую часть [2]:

3 Ay  u, Ay   a1 yx x  qy, x  , ( 6 )  1 y  x  0, x  , u  x   x 1  x / h12 2  x / h22 3  x / h32, x , ( 7 ) где где a1 h1, x2 , x3    0, x2 , x3 , x1  h1,  1  x  0, 2h1  x1  l1  2h1,  a1 l1, x2 , x3    l1, x2 , x3 , x1  l1  h1, a2  x1, h2 , x3     x1,0, x3 , x2  h2 ,   2  x  0, 2h2  x2  l2  2h2 ,  a2  x1,l2 , x2     x1,l2 , x3 , x2  l2  h2 , a3  x1, x2 , h3     x1, x2 , 0, x3  h3,  3  x  0, 2h3  x3  l3  2h3,  a3  x1, x2 , l3     x1, x2 , l3 , x3  l3  h3, x2 2 , x3 3; x1 1, x3 3; x1 1, x2 2.

 n1  D  R1  D1  D  R2 

 Ayn  f , Функция источников q, входящая в уравнение ( 6 ), носит  – образный характер и в подавляющем большинстве узлов сетки равна 0. Представим схему итерационного двухслойного модифицированного попеременно-треугольного метода в виде: yn1  yn ( 8 ) R2 y  31  ah1 yx  a2hx y  16 qy . (10) Операторы, введенные в соответствии с равенствами (9) и (10), являются сопряженными на множестве сеточных функций, обращающихся в ноль на границе сетки:

R1  R2 ,

A  R1  R2. Получены оценки для постоянных ∆ и δ:

  D  R1  R2 , R1D1R2   R1  R2 , .

4 а также сеточная функция d(x), определяющая элементы диагональной матрицы D. d  x  3  a1  h2 a2hx  q  , x  .

 1  6  h2 Будем считать, что v  x – решение краевой задачи a vx x  q v   a x  ,   1, 2,3; x  ,

q 3 2h 6 v  0, x  0, x  l . Положим c  x , x   max v  x, x   x1, x2 , x3  , x  c1  x2 , x3 ,  x2 , x3  2 3, если  1;  где c  x , x   c2  x1, x3 ,  x1, x3  1 3, если  2;  c3  x1, x2 ,  x1, x2  1 2 , если  3. Рассмотрим теперь другую совокупность трехточечных разностных задач:

a1 a wx x  q3 w   2 , w  0, x  0, x  l ,   1, 2,3,

h x  , x   x1, x2, x3 . Пусть b  x , x   max w  x, x   x1, x2 , x3   ,   1, 2,3,

x  b1  x2 , x3 ,  x2 , x3  2 3 , если  1; 2,3,если   1;   где b  x , x   b2  x1, x3 ,  x1, x3  1 3, если  2;  ,   1,3, если   2; b3  x1, x2 ,  x1, x2  1 2 , если  3; 1, 2, если   3.

 o 2   Тогда  d y ,1   A yo, yo ,   1.

    Положим q  max q  x1, x2 , x3 , x1,x2 ,x3

4  0x1l1h1  a1  x1, x2 , x3   a1  x1  h1, x2 , x3  ,  1  2 max a1 0, x2 , x3 , a1 l1, x2 , x3 , max

h1  2   2  42 max a2  x1, 0, x2 , a1  x1,l2 , x3 , max  a2  x1, x2 , x3   a2  x1, x2  h2 , x3  , h2  0x2 l2 h2  2  (9) (11) (12) (13) (14) (15) (16) (17)  3  42 max a3  x1, x2 , 0, a3  x1, x2 , l3 , max  a3  x1, x2 , x3   a3  x1, x2 , x3  h3  .

h3  0x3l3 h3  2  Введем постоянные величины (коэффициенты) в соответствии с равенствами: k1   /   q / 3,  1, 2,3. Выражение для определения функции d(x) запишется следующим образом:

1  d  x  31  ah21   b 2 a x  q6  b  c  kb1 2  .

 k1  2h

1 1   Оценка для параметра ∆ имеет вид:

  1 2   1,2,3 x  ,x   b  x , x  k1  x , x 2  c  x , x    ,   1, 2,3, x  , x  ,   max  max  (18) (19) 2,3, если   1,   ,   1,3, если   2,

1, 2, если   3. С учетом равенства δ=1, приходим к равенству, определяющему параметр 0 : 0  2 / . (20) В случае применения чебышевского ускорения итерационного процесса для числа итераций n0(ε), требуемого для достижения заданной точности ε справедлива оценка: 4   2    2   n0    ln  , n0    O  N0 ln   , N0  maxN1, N2 , N3 . (21) 2 2         4. Решение модельной задачи

Приведем постановку модельной задачи. Средняя проницаемость нефтеносных пластов составляет 10-12 м2. Это значение взято в качестве базовой проницаемости. Далее базовая проницаемость умножается на неоднородный коэффициент проницаемости, показанный на рис. 1 (сгенерирован случайно из диапазона 0.01…100 на поле 150×150 метров). Таким образом, разброс значений поля проницаемости составляет 4 порядка.

25 50 75 100 125

150 Рис.1. Распределение проницаемостей пласта Толщина пласта: 10 м., горизонтальные размеры области: 150×150 метров, граничные условия – непроницаемые стенки. Задавались суточные дебиты скважин, в одну из которых 100.00 51.79 26.83 13.89 7.20 3.73 1.93 1.00 0.52 0.27 0.14 0.07 0.04 0.02 0.01 (нагнетательную) вода закачивалась с дебитом 300 м3/сутки, а из другой (эксплуатационной) выкачивалась нефтеводоносная смесь с отрицательным дебитом 300 м3/сутки. Проницаемость пласта в зависимости от вертикальной координаты не менялась. Использовалась пространственно-трехмерная сетка с общим количеством шагов: 300×300×20, шаг по времени составлял 0.003 – 0,005 суток, сеточное уравнение для функции давления решалось усовершенствованным МПТМ. На рис. 2 приведены результаты численных экспериментов по расчету функции водонасыщенности.

Рис. 2. Распределение функции водонасыщенности через 10 (а), 20 (б), 50 (в), 200 (г) дней Результат расчета распределения функции давления через 50 дней приведен на рис. 3. В табл. 1. приведено количество итераций, требуемое для перехода с предыдущего слоя на новый слой по времени.

Таблица 1. Количество итераций, требуемое для перехода на новый слой по времени Номер временного слоя Результаты численных экспериментов демонстрируют значительное преимущество построенного варианта МПТМ по сравнению с употребительными методами методом Зейделя и методом верхней релаксации с шахматной упорядоченностью узлов – выигрыш в числе итераций составляет более чем 10 раз. 5. Параллельная реализация

Для параллельного решения модельной пространственно-трехмерной задачи фильтрации использовался установленный в Таганроге и эксплуатируемый с 2010 г. HPC-кластер Южного Федерального университета, который представляет собой многопроцессорную вычислительную систему с распределенной памятью, состоящую из 128 вычислительных узлов и общим числом ядер 2048. Опыт решения подобных сеточных задач на данной системе [3–8] показал, что декомпозиция по одному координатному направлению дает приемлемые значения эффективности при относительно небольшом числе, порядка 10-15, используемых вычислительных ядер. Поэтому была применена декомпозиция по двум координатным направлениям (рис. 4), при этом смежные области перекрывались двумя строками или столбцами узлов по направлению, перпендикулярному плоскости разбиения. Декомпозиция по трем координатным направлениям при существенной асимметрии в количестве узлов по горизонтальному и вертикальному направлениям не позволяет построить эффективный алгоритм. Рис 4. 2 D декомпозиция области Рассмотрим параллельный алгоритм расчета вектора поправки: (D  m R1)D1(D  m R2 )wˆ m1  rˆm , где R1 – нижне-треугольная матрица; R2 – верхне-треугольная матрица. Для вычисления вектора поправки нужно последовательно решить два уравнения: (D   R1) yˆm1  rˆm , (D m R2 )wˆm1  Dyˆm1 .

Вначале вычисляется вектор y m , при этом расчет начинается в левом нижнем углу. Затем из правого верхнего угла начинается вычисление вектора поправки wm . Схема расчета вектора y m изображена на рис. 5. На первом шаге вычислений первый процессор обрабатывает верхний слой. Затем осуществляется передача перекрывающихся элементов смежным процессорам. На следующем шаге первый процессор обрабатывает второй слой, а его соседи – первый. Передача элементов после расчета двух слоев первым процессором показана на рис.5.

Рис. 5. Схема обхода области для расчета вектора ym Результаты расчета ускорения и эффективности в зависимости от количества процессоров для параллельного варианта адаптивного попеременно-треугольного метода приведены в табл. 2.

Таблица 2. Зависимость ускорения и эффективности от количества процессоров Данный алгоритм позволил добиться ускорения более чем в 46 раз на 128 ядрах (по сравнению с последовательным алгоритмом), что в итоге позволяет в 300 – 2000 раз уменьшить время решения задачи по сравнению с традиционными последовательными алгоритмами – меКоличество процессоров 1 2 4 8 16 32 64 128 256 Ускорение

Эффективность 1 1,90 3,12 5,96 9,76 17,92 31,36 47,1 46,6 тодом Зейделя и методом верхней релаксации с шахматной упорядоченностью узлов при решении модельных задач с числом узлов до 3 200 000. 6. Заключение

Разработанный алгоритм позволил добиться ускорения более чем в 46 раз на 128 вычислительных ядрах (по сравнению с последовательным алгоритмом), что в итоге позволяет в 300 – 2000 раз уменьшить время решения задачи фильтрации по сравнению с традиционными последовательными алгоритмами – методом Зейделя и методом верхней релаксации с шахматной упорядоченностью узлов при решении модельных задач с числом узлов до 3 200 000, и в 10 – 60 раз по сравнению с их параллельными аналогами. Литература Parallel numerical implementation of spatially three-dimensional model of two-phase fluid filtration based on advanced alternating triangular method

A.I. Sukhinov1, A.E. Chistyakov2, A.A. Sukhinov3, L.A. Grigoryan4

Don state technical University1, Kalyaev Scientific Research Institute of Multiprocessor Computer Systems at Southern Federal University2, Skolkovo Institute of Science and Technology3, North-Caucasus Federal University4 The spatial-three-dimensional problem of two-phase filtration of incompressible fluid was being considered within the Buckley-Leveret model in view of the gravity and without capillary forces in variables "pressure-water saturation". The advanced modified alternating triangular iterative method (MATM) that takes into account the availability of sources (sinks) in a relatively small number of nodes with delta-shaped (generalized) functions was considered for the discrete model that approximates the boundary-initial filtration problem.

The numerical solution of computationally laborious task in the layer with significant heterogeneity of permeabilitywas performed, which are changed to four magnitude orders within the layer with a significant advantage in required iterations to achieve a given accuracy − in 10 and more times on grids with up to 3.2 million nodes. The parallel algorithm of numerical implementation of the developed sequential variant of advanced MATM was proposed in this paper which are based on the decomposition method in two spatial directions. It is performed on a multiprocessor system with shared memory − the computer system of the Southern Federal University in Taganrog that are contained 2048 cores with the peak performance 18.8 teraflops. The significant advantage of the proposed parallel algorithm was shown, compared to traditionally-used methods for numerical solution of these problems − the time simulation of the flooding process of layer system was reduced in tens of times.

Keywords: filtration problem, water saturation, modified alternating triangular method.

References

1. Chetverushkin B. N. , Morozov D. N. , Trapeznikova M. A. , Churbanova N. G. , Shil' nikov E. V. An explicit scheme for the solution of the filtration problems // Mathematical Models and Computer Simulations . December, 2010 . Vol. 6 , No. 2. P. 669 - 677 .

2. Sukhinov A.I. Modifitsirovannyy poperemenno - treugol'nyy metod dlya zadach teploprovodnosti i fil'tratsii [Modified alternating triangular method for heat-conduction and filtration problems] // Rostov-na- Donu , RSU, Vychislitel'nyye sistemy i algoritmy [Rostov-on- Don , Rostov State University, Computer systems and algorithms] . 1984 . P. 52 - 59 .

3. Sukhinov A.I. , Chistyakov A.E. Parallel 'naja realizacija trehmernoj modeli gidrodinamiki melkovodnyh vodoemov na supervychislitel'noj sisteme [Parallel implementation of threedimensional hydrodinamics model of shallow waters on a supercomputer system] . Vychislitel'nyye metody i programmirovaniye [Numerical Methods and Programming] . 2012 . Vol. 13 . P. 290 - 297 .

4. Sukhinov A.I. , Protsenko E.A. , Chistyakov A.E. , Shreter S.A. Sravneniye vychislitel'nykh effektivnostey yavnoy i neyavnoy skhem dlya zadachi transporta nanosov v pribrezhnykh vodnykh sistemakh [Comparison of computational efficiencies of explicit and implicit schemes for the sediment transport problem in coastal zones] . Vychislitel'nyye metody i programmirovaniye [Numerical Methods and Programming] . 2015 . Vol. 16 , No. 3. P. 328 - 338 .

5. Sukhinov A.I. , Chistyakov A.E. , Semenyakina A.A. , Nikitina A.V. Parallel 'naya realizatsiya zadach transporta veshchestv i vosstanovleniya donnoy poverkhnosti na osnove skhem povyshennogo poryadka tochnosti [Parallel implementation of transport problem of substances and restoration problem of the bottom surface on the basis on schemes of high-order of accuracy] . Vychislitel'nyye metody i programmirovaniye [Numerical Methods and Programming] . 2015 . Vol. 16 , No. 2. P. 256 - 267 .

6. Sukhinov A.I. , Chistyakov A.E. , Protsenko E.A. Matematicheskoye modelirovaniye transporta nanosov v pribrezhnykh vodnykh sistemakh na mnogoprotsessornoy vychislitel'noy sisteme [Mathematical modeling of sediment transport in coastal water systems on multiprocessor computer system] . Vychislitel'nyye metody i programmirovaniye [Numerical Methods and Programming] . 2014 . Vol. 15 . P. 610 - 620 .

7. Sukhinov A.I. , Nikitina A.V. , Chistyakov A.E. , Semenov I.S. Matematicheskoye modelirovaniye usloviy formirovaniya zamorov v melkovodnykh vodoyemakh na mnogoprotsessornoy vychislitel'noy sisteme [Mathematical modeling of the formation of suffocation conditions in shallow basins using multiprocessor computing systems] . Vychislitel'nyye metody i programmirovaniye [Numerical Methods and Programming] . 2013 . Vol. 14 . P. 103 - 112 .

8. Sukhinov A. I. , Nikitina A.V. , Chistyakov A. E. Numerical simulation of biological remediation Azov Sea / / Mathematical Models and Computer Simulations . 2012 . Vol. 24 , No. 9. P. 3- 21 .