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

УДК 519.622.2

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ФИЛЬТРАЦИИ НА МНОГОПРОЦЕССОРНЫХ СИСТЕМАХ

М, В, Васильева

В работе рассматривается изотермическая однофазная одноком-понентная фильтрация флюида. Уравнение фильтрации получается комбинированием закона сохранения массы и закона Дарси, которые детально обсуждаются во многих книгах [1,2]. Полученное нестационарное нелинейное дифференциальное уравнение в частных производных с граничными и начальными условиями аппроксимируется стандартными конечно-разностными аналогами по пространственным переменным, что приводит к системе дифференциальных и алгебраических уравнений.

1. Математическая модель

Для моделирования изотермической фильтрации флюида в трехмерном пласте [1] запишем закон сохранения массы для трехмерного течения в пористой среде:

I. {рф) = -У(ри) + др, (1)

где р — плотность флюида, ф — пористость, и — скорость фильтрации и д — отбор в единицу времени. Добавим к нему соотношение между скоростью фильтрации и градиентом давления. Дифференциальная форма этого соотношения имеет вид

к

и=— (УР + РЕ), (2)

© 2010 Васильева М. В.

где к — тензор абсолютной проницаемости пористой среды; р — вязкость флюида; g — вектор ускорения свободного падения.

Уравнение однофазной фильтрации получается комбинированием закона Дарси и уравнения сохранения массы. Плотность флюида является функцией от давления:

др.

(3)

Для получения более удобной формы предположим, что сжимаемость флюида и порового пространства, определяемые выражениями

1 ЗУ

1 Эр р др

1 дф

г ф др

(4)

постоянны во всем интересующем нас диапазоне давлении.

Выражая член с производной по времени через др/дЬ, получим д др др дф др др т{рф) = ФТРЪ + = Мч + Сг) д?

(5)

(6)

др рк рр(сг + сг)— = V (— (Ур + р%) I + др.

Дополним уравнение (6) условием непротекания на границе:

ип = О,

и заданным начальным распределением:

р{х,у, г,Ь = 0)= р0.

Таким образом, мы получили начально-краевую задачу нестационарной нелинейной фильтрации, описываемую нелинейным дифференциальным уравнением в частных производных с переменными коэффициентами:

С(р)^=У(\(р)(Чр + р%)) + др, (7)

где

С(р) = фр(с/ + Сг), \(р) = Поставленную задачу будем решать численно.

рк р '

2. Дискретизация по пространственным переменным

Для численного решения задачи (6) проведем стандартную конечно-разностную аппроксимацию пространственных переменных на равномерной сетке [3] и получим систему обыкновенных дифференциальных уравнений (ОДУ):

f-i dpi j fc __PiJjk Pi,j,k Pi—l,j,k

\ i^jJ+ljfc -» Pi,j,k Pi,j — l,k

+ \,j+1/2,к д^ \,j-l/2,k д^

к - Pi,j,k-1

+ \,j,k+1/2-д"^--\j,k-l/2-д^з-

pi,j,k+l pi,j,k л pi,j,k Pi,j,k-1

+ \,j,k+1/29-д^--\,j,k-l/29-д^-

+ A xly^ (®)

где q = qAi,ji,k и 5i}j}k — символ Кронекера, который всюду равен О, кроме точечных источников, где он равен 1.

Добавим к уравнению (8) дискретные аналоги граничных условий и получим систему дифференциальных и алгебраических уравнений (DAE, diffirencial and algebraic equation), которую можно записать в общем виде:

p(0)jk= P>. (10)

Специфику рассматриваемой задаче придает наличие скважин (радиуса rw), диаметр которых неизмеримо мал по сравнению с размерами пласта. Для учета скважин в численных моделях фильтрации основываемся на допущении того, что вблизи скважины течение описывается аналитическим решением. Этот подход впервые предложен Писма-ном [4]. и развит в работах А. Н. Чекалина [5] и Д. В. Шевченко [6]. Поскольку скважины очень малы, рассмотрим их как линейные стоки/источники для скважин, полностью вертикально пронизывающих

пласт. Предположим что вблизи скважин характер течения близок к радиальному. Запишем радиальный закон Дарси:

Q = 2nrAz^. (11)

¡л dr

Перепишем его в следующем виде и проинтегрируем:

rw rw

dr dr r

dp = WI—, dp = WI —, Ap(r) = Wllnl — r J J r \rw

rr

p(r) =Pwf -Wlln^J-^j,

где

wi= "g

(12)

2nkAz'

Запишем конечно-разностную аппроксимацию уравнения материального баланса:

~^{Рг+1,3,к ~ Pi,j,k) - (P¿,j,fc -Pi-l,j,k) ~(Pi,j+l,k — Pi,j,k) — ~T 9 (Pi>,j,k — Pi,j-l,k)

Q =0. (13)

AxAyAz

уу

Для выполнения уравнения (13) в узлах вблизи скважины введем в него поправочные коэффициенты и вычислим их подстановкой давлений, определяемых уравнением (12). Получим поправочные коэффициенты в половинных узлах:

п

Г1гд±1/2,]д,кдТГ1гд,]д±1/2,кд = Т^^Д^'

1

'Ъд±1Лд±1/2,кд,'Шд±1/2^д±1,кд -

1

ГЦд±1^д±3/2,кд,ГЦд±3/2^д±1,кд ~ ^^ 5 '

7Г -2

ЦЪд±Ъ12,Зд,кд1ЦЪд,Зд±Ъ12,кд ~

которые вводятся в шаблон для узлов (гч ± kq), (гч± 1, kq), (гч ± ± 1, kq) вблизи скважины (iq,jq, А^ = Ау.

3. Методы расчета

Для численного решения DAE на многопроцессорной системе воспользуемся пакетом IDA из библиотеки параллельных алгоритмов SUNDIALS [7]. Код IDA это реализация на С фортрановской библиотеки DASPK. IDA решает задачу Коши для системы дифференциальных и алгебраических систем, заданных в общем виде:

4»|)=о. ,№■) = *, =

Для его решения применяется метод дифференцирования назад (BDF, backward differencial formulae) [8]:

= Ь-п-^-, (14)

¿=0

где Нп = Ьп — ¿п_1 и ап^ известны. После применения метода ВОГ к системе ОАЕ получается система нелинейных алгебраических уравнений, которая решается на каждом временном шаге:

О{уп) = ^(^п, Уп, К1 ^(ап^= 0- (15)

Для решения (15) применяется метод Ньютона, который приводит к системе линейных уравнений для каждой итерации Ньютона, в форме

Л (УТ1 — У?) = —0(у?), (16)

где Л = ¿О/ёу — якобиан.

Полученная система линейных алгебраических уравнений решается методами подпространства Крылова с предобуславливанием [9].

4. Результаты расчетов

Расчеты проводились в области

[0, 2000 м] х [0,2000 м] х [—2500 м, —2000 м].

stress 2850.6

Рис. Распределение давления при откачке (с^) 500 и 1000 баррелей нефти в день. Поверхность для давления 2500 ней.

6250 6240

6230

stress 6285.3

6259.1

6232.8

6206.6

6180.3

Рис. 2. Распределение давления при откачке (с^) 500 и 1000 баррелей нефти в день. Срез но г=-2220 м.

stress 2654.9

2563.7

2472.6

12381.5

2290.3

Рис. 3. Распределение давления для двух нагнетательных скважин (q = 1000 баррелей воды в день). Срез ио z=-'2'2'20 м.

При этом скважины полностью пронизывали пласт и располагались в ячейке x = МО м, y = МО м и ж = 1500 м, y = 1500 м. Расчеты проводились на кластере СКЦ ИМИ ЯГУ при разных количествах процессов для различных сеток. Расчитывалось распределение пластовых давлений для различных флюидов (вода, нефть) на 10 лет. Результаты расчетов проиллюстрированы на рис. 1 3.

При этом использовались следующие параметры модели: kx = ky = kz = 100 мДарси, ф = 0.2, cr = 6.0 * 1О-6 1/пси, pwater = 1000, кг/м3, cwater = 3.0 * 1О-6 1/пси, роп = 611.7, кг/м3, с0п = 22.14 * 10-6 1/пси.

ЛИТЕРАТУРА

1. Aziz К., Sett ату A. Petrolium reservoir simulation. Amstertam: Elsevier Appl. Sci. Publ., 1979.

2. Reservoir simulation manual. Heriot Watt University. Department of Petroleum Engeneering, 2002.

3. Самарский А. А., Николаев E. С. Методы решения сеточных уравнений. М.: Мир, 1983.

4. Peaceman D. W. Interpretation of well-block pressures in numerical reservoir simulation. SPERE 5. 1990.

5. Чек алии A. H. Численные решения задач фильтрации в водонефтяных пластах. Казань: Изд-во КРУ, 1982.

6. Шевченко Д. В. Применение многосеточных методов для расчета давления в нефтяном пласте // Мат. моделирование. 2002. Т. 14, № 8. С. 113^118.

7. Sundials - www.llnl.gov/casc/sundials

8. Hairer Е., Noersett G., Wanner S. P. Solving ordinary differential equations. V. 1. NonstifF problems, 2 ed. Springer, 2008.

9. Saad Y. Iterative methods for sparse linear systems. 2 ed. Cambridge: Cambridge Univ. Press, 2003.

г. Якутск

31 мая 2010 г.