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

УДК 519.63

ПАРАЛЛЕЛЬНОЕ ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ЗАВОДНЕНИЯ НЕФТЯНОГО МЕСТОРОЖДЕНИЯ*)

II. \I. Афанасьева, М. В, Васильева, П. Е, Захаров

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

Модель представляет собой начально-краевую задачу для определения полей насыщенности, давления и температуры. Для численного решения проводится стандартная конечно-разностная аппроксимация, которая приводит к разреженной системе нелинейных уравнений для каждой неизвестной. Численное решение таких систем является достаточно ресурсоемкой задачей и требует применения вычислительных кластеров. Для параллельной реализации задачи в работе применяется свободно распространяемый пакет PETSc (Portable Extensible Toolkit for Scientific Computing) [5]. Он включает в себя средства для

*) Работа выполнена при поддержке фцп «Научные и научно-педагогические кадры инновационной России» на 2009-2013 г.г. (гос.контракт № 02.740.11.0041)

@ 2011 Афанасьева Н. М., Васильева М. В., Захаров П. Е.

решения систем иелинеиных и линеиных уравнении с использованием метода Ньютона и методов подпространств Крылова с предобуслов-лпванпем [6]. Для организации параллельной работы алгоритмов на многопроцессорных системах с распределенной памятью в пакете применяется метод декомпозиции расчетной области, а в качестве коммуникационной библиотеки используется интерфейс передачи данных MPI (Message Passing Interface) [7].

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

Рассмотрим математическую модель неизотермической двухфазной фильтрации в однородном пласте. Предполагаем, что одна из фаз больше смачивает пористую среду, чем другая. Индексы w и n обозначают соответственно смачивающую и несмачиваюгцую фазы. Запишем уравнения неразрывности для каждой из фаз:

д

■7^(<ppisi) + V(piUi) = piqi, l = w,n, (1)

где pi — плотность l-й фазы, ф — пористость, щ — фазовая скорость фильтрации, qi — объемные источники/стоки и si — пасыщенность 1-й фазой.

Предполагаем, что двухфазная система занимает все поровое пространство:

s^ sn = 1, (2)

каждая фаза движется в пористой среде по обобщенному закону Дарси

k

щ =--l = w,n, (3)

Pi

где k — абсолютная проницаемость среды, pl — коэффициент динамической вязкости l-й фазы и krl — относительные фазовые проницаемости, характеризующие зависимость скорости фильтрации от насыщенности, а Ф| — потенциал l-й фазы, который выражается следующей формулой:

= pi + pi gz. (4)

Здесь ось ^ направлена вертикально вверх, р1 — давление /-й фазы, д — ускорение свободного падения.

Из-за кривизны и поверхностного натяжения между границами раздела двух фаз давление смачивающей фазы меньше, чем несмачи-вающей. Эта разность давлений определяется как капиллярное давление:

Рс= Рп - Рш- (5)

Эмпирически капиллярное давление является функцией от насыщенности.

Введем \ как фазовую подвижность:

АI = —, I = ги,п, (6)

и Хт как суммарную подвижность:

Хт = Хт + Хп- (7)

Определим также функцию

/г = "Г", 1 = ш,п, (8)

Хт

и суммарную скорость

ит = иш+11п. (9)

Сложим уравнения сохранения массы для смачивающей и несмачива-ющей фаз и получим следующее уравнение:

Уит = + Яп = Ят - (Ю)

Выразим иш через скорость песмачивающей фазы:

= + кХтУрс + к(рп - рги)Хгиде3, (11)

Хп

где е3 = (0,0,1)т.

Поскольку иш = ит — ип, то ип выразим через суммарную скорость:

Цп = /п^т — к/тХп( Урс + {рш — Рп) дез)- (12)

Подставим (12) в уравнение неразрывности для несмачивающей фазы:

dsw

ф-7^- + v (fwuT + kfwXn(Vpc + (рп - pw)ge3)) = qw. (13)

Преобразуем (13) и получим уравнение для определения насыщенности:

Ф-^jf + UTfLVs* + V[kfwX„p'cVsw} + к(рп - pw)gd= qw - fwqT-

Z (14)

Так как рассматриваемый процесс неизотермический, т. е. температура в пласте не остается неизменной, запишем систему уравнений, описывающую теплопередачу при двухфазном течении флюидов в пористой среде [8,9].

Уравнение теплового баланса для пористой среды (solid):

dT

csps—^ = V(ksVTs). (15)

dt

Уравнение теплового баланса для воды:

(dT \

cwPw (+ V(uwTw)j = V(kwVTw). (16)

Уравнение теплового баланса для нефти:

СпРп + V(unTn)^j = V(knVTn). (17)

Предположим, что фазы находятся в локальном термодинамическом равновесии, значит, температуры всех фаз и скелета пористой среды в каждой точке пласта совпадают, Tw = Tn = Ts. Тогда, сложив уравнения теплового баланса, получим уравнение теплопроводности в пористой среде с учетом фильтрации двухфазной жидкости:

dT

[(1 - 4>)csps + 4>{swCwpw + (1 - Sw)cnpn)} —

[sw cw Pw^-w ""b (1 sw) cnPnHn] VT

= V([(l - ф)кв+ ф{Swkw + ( 1 - Sw)kn)}VT) + q'", (18)

где Т — температура, сш,сп,с8 — удельные теплоемкости, кш,кп,к8 — теплопроводности флюидов и пористой среды соответственно, я'" — удельная мощность источников тепла.

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

Уит = ят, (19)

дЗш

Ф+ и т/^У««, + У^/иД«^ Ув«,]

+ к(р„ - = Яги - /ш<7т, (20)

дТ

[(1 - Ф)с8Рз + Ф^шСшРш + (1 - ««>)с„/Э„)] —

+ [Зш сшршИш ""Ь (1 Зш)спрпИп]УТ

= У([(1 — ф)кв+ ф{8ш кш + ( 1 — Зш) кп)] УТ) + Я''' , (21)

Для этой системы задаются начальные условия, а также на внешней границе ставятся условия непротекания.

Будем предполагать, что через нагнетательные скважины закачивается однородная подогретая жидкость с заданным дебитом и температурой Ть. Входной и выходной поток равны сумме объемных источников и стоков, ят = Яш + Яп■ На вход подается только вода и яш = ят, Яп = 0. На выходе первое время вытекает нефть, а затем смесь воды и нефти, яш = ятХш/Хт, яп = ЯтХп/Хт [2,4]. Фазовые потенциалы Ф; определяются в соответствии с формулировкой через фазовое давление, р = рп:

УФ ш = Ур — р'с Узш+ Ршде3, УФ п = Ур + Рп^з-

Для определения капиллярного давления, функциональной зависимости между относительными проницаемостями фаз и насыщенностью применяем модель Вгоокв—Согеу [10]:

7 2 + 3в 1 /Л \2 /-1 2 + 3в •

кги1 = ,зе » , кгп = (1 — ве) (1-ве о ),

где зе — эффективная насыщенность,

— "I $гги ^:: $11} Я

1 Згш Згп

ив — постоянный коэффициент. Капиллярное давление зависит от насыщенности и входного давления ре:

Рс(йп,) =

Коэффициенты динамической вязкости от температуры будем вычислять по формуле Булыгина для нефти и формуле Пуазейля для воды [11]:

Рп(Т) = 0-001/ (0^2 + 0-005 (Т — 273-15))2 ,

Рш(Т) = мо/(1 + 0-0337(Т — 273-15) + 0-000221(Т — 273^)2), где мо — вязкость воды при температуре 273.15 К.

2. Дискретный аналог и алгоритм решения

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

Рассмотрим задачу (19)-(21) в параллелепипеде со сторонами Ьх, Ьу, Ьг. Введем равномерную сетку с шагами Нх,Ну,Нг по направлениям х, у, г и координату то времени разобьем па временные шаги т. Положим з = зш. Производная по времени Аз^/АЬ аппроксимируется обратной разностью

а \ п+! (Зп+1 _

ааг,],к \ ^ У г,],к "г,],к

Л

37^1 - зп.

(22)

Запишем конечно-разностную аппроксимацию уравнения для насыщенности (20):

Ф^ + и/^Ув + У[/г/шКр'сУй\ + к(рп - = - (23)

(е^1 - оП )

фУК^кт 1'3'к) + (В1+1/2ГЫ-

+ (В— /2)+ {ег,з,к - в— ¿,к)п+1 + (в^!^) — в1,],к)п+1 + (В— /ъ)+ — $1,3-1 ,к)п+1

+ вк+1/2) — $13,к)П+1 + В— /2) + {$1,3,к — $1,з,к-1 )П+1

+ (kXnfwP' c)i+l/2 {$1+1,з,к — $1,з,к)П+1 Ну Нг

— (k\nfwp' с)—1 /2 — в— ¿,к)п+1 Ну Нг + {k\nfwP'c)j+l/2 (Ч3+1,к — КНг

— (кХп fwр'с)з —1 /2 (si,j,k — si,j-l ,к)п+1 Нх Нг + (kXnfwP'с)к+1/2 (si,j,^+l — '¡¿3,к)п+1 НхНу

— (к\п fwP'c)k-l/2 (si,j,k — )П+1 ЬхНу

+ {к{рп — рw) З) i,j,k( fw^nk+l/2 — fwХnk — l/2)П^ НхНу

= — fw i,j,k , (24)

где

Г о, в > о,

В—

1 В, В < о, ^ ;

Г В, В > о,

В

I о, В < о,

■#¿+1/2 = и^+1/2,3,к( fwi+l,j,k — fwilj,k) / — si,j,k) Ну Нг, (27)

В— /2 = UTi— /23,к{ fwi,j,k — fwi-l 3,к)/(si,j,k — З,к) Ну Нг . (28)

Аналогичным образом определяются Вз+1/2, Вз—1 /2, Вк+1/2, Вк—/2. Аппроксимация функции Хт в промежуточных пространственных координатах происходит по схеме со «взвешиванием вверх по потоку» [2].

Дискретный аналог уравнения суммарной скорости (19) запишется следующим образом:

(и^+1/2 — — 1 /2) Ну Нг + (итз+1/2 — ит— /2) Нх Нх

+ (итт/2 — -От—/2) Ну Нх = . (29)

Уравнению температуры (21) поставим в соответствие неявную разностную схему:

Тп — Тп

[(1 - ф)с3р3 + ф(вси]ри] + (1 - 3)спрп)]У-^-—-

т

+ Ьг+1/2 — Т1,з,к)п+1 + Ьг—1/2 {Т1,],к — Тг—1 ,],к)п+1

+ ^¿+1/2 {Т1,]+1,к — Т1,],к)п+1 + /2 {Тг,3,к — Тг,]~1 кУ^ + Ьк+1/2 {Т1,],к+1 — Т1,],к)п+1 + Ьк — 1 ¡2 {Тг,3,к — Тг,],к — 1 )п+1

+ Аг+1/2р'ы+1/2 (Зг+1,],к — Зг^,к) + ¿1/2р^+1/2 (Зг,]+1,к — Зг^,к) + Ак+1/2 (р'ок+1/2 (3г,3,к+1 — Зг,],к) — (Рп — Я)

= кг+1/2 — Т1,],к)п+1 Ьу Нг — к—1 /2 {Т1,],к — Тг—1 ,],к)п+1 Ьу Нг

+ ¿¿+1/2 {Т1,]+1,к — Т1,],к)п+1 ЬхЬг — к] — 1/2 {Тг,3,к — Тг,]—1 ,к)п+1 ЬхЬг + кк+1/2 {Т1,],к+1 —Т1,],к)п+1 ЬхЬу — кк — 1 /2 {Тг,],к —Тг,],к — 1 )п+1 НхЬу + ({"^ ^ к ,

(30)

где

Ът±1 /2 = ((зСЫРwfw + (1 — з)СпРп1п)ит)±/2 , т = 1,3,к Ат+1/2 = ((зcwРw — (1 — з)спРп)kfw^п)п^1/2, т = ¿,3, к

кт±1 /2 = ((1 — Ф)кв+ ^зк^ + (1 — з)кп))п±1/2, т = 1,3, к.

Для численного решения поставленной дискретной задачи получаем следующий процесс. На каждом временном слое рассчитываем давление при фиксированной насыщенности и температуре (значения берутся с предыдущего временного слоя). Затем определяем насыщенность на текущем слое с учетом вновь найденных значений давления. Далее вычисляем температуру, после чего осуществляем переход на следующий временной слой. Заметим, что при аппроксимации по времени насыщенности могут определяться как явно, так и неявно, а давления и температуры вычисляются по неявной схеме [12].

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

пространяемый пакет PETSc (Portable Extensible Toolkit for Scientific Computing).

PETSc содержит набор средств для параллельного (а также и последовательного) решения дифференциальных уравнений с частными производными. Эти средства включают в себя решатели линейных и нелинейных уравнений, объекты для управления распределенными данными (DA) и организацией обмена между ними, где в качестве коммуникационной библиотеки используется интерфейс передачи данных MPI (Message Passing Interface).

Решение нелинейных систем происходит с использованием модуля SNES (Nonlinear Solvers) для уравнений вида

F(u) = 0, (31)

которые решаются методом Ньютона со следующим итерационным процессом:

м„+1 =un + Sun. (32)

Здесь 5un находится из решения системы линейных уравнений:

J(uri) 5un = -F(uri), (33)

где J — якобиан функции F. Поскольку вычисление якобиана является сложной задачей, его вычисляют, используя следующие дифференциальные разности:

Jj Ю = AUj F(un) « [F{un + hej) - F(un)]/h, (34)

Для решения системы линейных уравнений используется семейство итерационных методов подпространств Крылова (KSP: Krylov subspace method) с предобуславливанием (PC: preconditioner).

3. Результаты

Приведем результаты параллельных численных расчетов, проведенных на вычислительном кластере ИМИ СВФУ. Кластер состоит

Рис. 1. Распределения давления. МПа.

Рис. '2. Распределения водоиасыщешюсти.

Рис. 3. Распределения температуры, К.

из четырех вычислительных узлов, в каждом узле имеются два че-тырехъядериых процессора Intel Xeon Е5450 (3.00 GHz) с оперативной памятью 16 Gb.

Расчеты проводились со следующими входными данными: Lx = 1000м, Ly = W00м, Lz = 50м, k = 1 0~13м2, ф= 0.2, ^ = 1ЛПа-с, pw = W00кг/м3, рп = 900кг/м3, ps = 2500кг/м3, kw = 0.6Вт/(м-К), kn = 0.15Вт/(м-К), ks = 1.ШВ^мCw = 4187Дж^г^, Cn = 2430Дж/(кг-К), c^00 Дж^г^^^ка 100 * 100 * 20, количество запущенных параллельных процессов 16. Относительные проницаемости определяются по модели Брукс — Корей с параметрами srw = 0.2, srn = 0.15, в = 2, pe = 104. В начальный момент времени t = 0 давление, водонасыщеность, температура распределены равномерно: po = ЮМПа, so = 0.2, Tq = 290 К. Для решения системы линейных уравнений применили метод GMRES с мультигрид предобуславлива-телем.

На рис. 1 3 представлены распределения давления, водонасыщен-

t

да. Четыре добывающие скважины расположены по бокам в верхней части области, с одинаковыми дебитами, д = 21.6 м3/сутки, и одна нагнетающая скважина посередине в нижней части с дебитом д = 86,4 м3/сутки.

Рис. 4,5 иллюстрируют динамику фронта давления и водонасы-хг

скважина, справа одна добывающая скважина. Здесь наблюдается прорыв воды к добывающей скважине.

На рис. 6 показан график зависимости от времени дебита извлекаемой нефти дп.

Таким образом, с использованием пакета РЕТБс численно реализована трехмерная задача неизотермической двухфазной фильтрации с учетом капиллярного давления.

0@ 3 месяца

ее 3 года

Рис. о. Распределения водонасыщешшсти при прорыве воды.

100 200 300 400 500 600

Бремя, сутки

Рис. 6. Добиты добывающей скважины.

ЛИТЕРАТУРА

1. Жолтов Ю. 11. Разработка нефтяных месторождений. М.: Недра. 1998.

2. Азия X.. Сеттарм 9. Математическое моделирование пластовых систем. М.; Ижевск: Институт компьютерных исследований. '2004.

3. Chen Z.. Guanrcn Н.. Yuanlc М. Computational methods for multiphase Hows in porous media. Society for Industrial and Applied Mathematic. 2006.

4. Коновалов А. Н. Задачи фильтрации многофазной несжимаемой жидкости. Новосибирск: Наука, 1988.

5. Balav S., Buscbelman К., Gropp W. D., Kausbik D., Kneplev M. G., Mclnnes L. C., Smith B. F., Zhang H. PETSc Web page, http://www.mcs.anl.gov/petsc

6. Saad Y. Iterative Methods For Sparse Linear Systems. Sec. Ed. Cambridge: Cambr. Univ. Press, 2003.

7. MPl: A message-passing interface standard. International J. Supercomputing Applications. 1994. V. 8, N 3/4.

8. Эккерт Э. P., Дрейк P. M. Теория тепло- и массообмена. М.-Л.: Госэнергоиздат, 1961.

9. Лыков А. В. Теория теплопроводности. М.: Высш. школа, 1967.

10. Brooks R., Corey A. Hydraulic properties of porous media. Vol. 3 of Colorado State University Hydrology Paper. Colorado State University, 1964.

11. Васильев В. И., Попов В. В., Тимофеева Т. С. Вычислительные методы в разработке месторождений нефти и газа. Новосибирск: Изд-во СО РАН, 2000.

12. Вабигцевич П. Н. Явно-неявные вычислительные алгоритмы для задач многофазной фильтрации // Мат. моделирование. 2010. Т. 22, № 4, С. 118^128.

г. Якутск

14 февраля 2011 г.