Software & Systems

Основным методом разработки нефтяных месторождений в России является заводнение. Большинство нефтяных месторождений находятся на третьей или четвертой стадии разработки. Следовательно, в продукции добывающих скважин присутствует большое количество воды. В этих условиях инженерам по разработке месторождений необходимо решать много задач по контролю и регулированию разработки. Для быстрого оценочного решения подобных задач может быть использована модель на базе концепции суперэлементов. Из-за большой неоднозначности исходных данных необходимо адаптировать модель к фактически замеренным показателям разработки (дебит нефти, газа и воды; забойное давление; пластовое давление и пр.). Задача адаптации модели сводится к решению задачи оптимизации. Кроме того, с помощью теории оптимизации и математического моделирования решают задачи оптимизации разработки нефтяных месторождений (оптимизации заводнения). Данная работа посвящена подходу, позволя- ющему с помощью простых методов (метод Нью- тона, метод сопряженных градиентов) решать задачи оптимизации при суперэлементном моделировании разработки нефтяных месторождений. Используя только Excel Microsoft, описываемым методом можно решать задачи оптимизации с количеством параметров до 16. Используя специальные программы, количество параметров можно увеличить.

Математическая модель двухфазной фильтрации нефти и воды

Система уравнений, описывающая упругую двухфазную фильтрацию, состоит из двух уравнений сохранения количества основных компонентов (объемов нефти и воды в стандартных условиях) и обобщенного закона Дарси, а капиллярными и гравитационными силами пренебрегается [1–6]:

– уравнения сохранения объемов нефти и воды

,                             (1)

;                    (2)

– обобщенный закон Дарси

,                                        (3)

,                                       (4)

где m –  пористость; So и Sw, Bo и Bw,  и ,  и , kro и krw, mo и mw – насыщенность, объемный коэффициент, вектор скорости фильтрации, плотность источника (стока), моделирующего работу скважины, относительная фазовая проницаемость, вязкость нефти и воды соответственно;  – плотность источника (стока), моделирующего работу законтурной области; k – абсолютная проницаемость; P –  пластовое давление.

Систему уравнений (1)–(4) необходимо дополнить замыкающими соотношениями:

So + Sw =1,                                                           (5)

m = m0[1 + Cr(P – P0)],                                      (6)

Bo = Bo0[1 – Co(P – P0)],                                          (7)

Bw = Bw0[1 – Cw(P – P0)],                                   (8)

где m0 – пористость при начальном пластовом давлении; Bo0 и Bw0 – объемный коэффициент нефти и воды при начальном пластовом давлении; Cr, Co и Cw – сжимаемость пласта, нефти и воды; P0 – начальное пластовое давление.

Систему уравнений (1)–(4) с учетом соотношений (5)–(8) можно преобразовать к следующей системе дифференциальных уравнений для водонасыщенности S º Sw и давления P [7]:

                       (9)

                           (10)

Уравнения (9) и (10) дополняются начальными условиями P = P(x, y, z, t = 0), S = S(x, y, z, t = 0) и граничными условиями непротекания на внешних границах. Для моделирования законтурной водоносной области используются источниковые слагаемые  в уравнениях (9) и (10) [1].

Численная схема

Рассмотрим расчетную сетку, составленную из так называемых суперэлементов [8, 9] – ячеек с размерами, в плане сопоставимыми с расстоянием между скважинами (300–500 м). Тогда количество ячеек в модели будет примерно равно количеству скважин. Скважины используются как центры ячеек. В неразбуренной части объекта могут быть введены фиктивные скважины для построения сетки. Использование такой сетки позволяет сократить вычислительные затраты в тысячи раз [8].

Задача построения суперэлементной сетки в плане (2D-сетки) сводится к задаче построения для каждой скважины зоны дренирования, которая бы­ла введена А.П. Крыловым, а она, в свою очередь, к построению областей Вороного [7] (рис. 1):

-     нанесение на карту скважин и границ;

-     построение выпуклой оболочки по внешним скважинам;

-     построение триангуляции Делоне по скважинам;

-     построение областей Вороного по триангуляции.

Для преобразования суперэлементной сетки в 3D-сетку построенную в плане сетку копируют для разных слоев, изменяя при этом глубину залегания кровли ячеек и их толщину.

Большой размер ячеек позволяет использовать полностью явную схему [10]. Рассмотрим такую схему для уравнений (9) и (10) на сетке суперэлементов в плане (2D-плоская модель):

где n – номер шага по времени; i – номер ячейки; Dtn+1 – шаг по времени; Fi, hi и (m0)i – площадь, толщина и пористость при давлении P0 ячейки i;  – дебит фазы a = o, w скважины в ячейке i на шаге по времени n;  – приток воды из законтурной области в ячейку i;  – переток фазы a = o, w между ячейками i и j на шаге по време- ни n; Ni – количество соседних ячеек для ячейки i.

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

Переток между ячейками i и j (рис. 2) определяется как , a = o, w, где Aij – длина совместного ребра ячеек i и j; Lij – расстояние между центрами ячеек (скважинами) i и j; ;

; .

Моделирование скважин

Дебиты скважин по жидкости ql, нефти qo, воде qw, обводненность W и забойное давление Pw опре- деляются из следующих соотношений:

,

,

,

,

где  – коэффициент продуктивности скважины i по жидкости, который определяется как

,

где (rw)i – радиус скважины по долоту;  – скин-фактор скважины; (Rc)i – эквивалентный радиус блока, определяемый как .

Для добывающих скважин необходимо задать дебит жидкости или забойное давление, для нагнетательных – расход воды или забойное давление.

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

где (Sc)i – водонасыщенность на фронте вытес- нения; (Swcr)i – насыщенность связанной водой; (Sowcr)i – насыщенность остаточной нефтью; коэффициенты Ai, Bi, ai, bi определяются по результатам интерпретации исследований керна и могут зависеть от пористости или проницаемости.

Водонасыщенность на фронте вытеснения (Sc)i  определяется из соотношения

,

где f(Si) – функция Бакли–Леверетта:

.

Предложенная 2D-расчетная схема легко обобщается на случай 3D. При этом приток в скважину из каждого слоя определяется согласно его продуктивности при известном дебите скважины по жидкости или забойном давлении.

Решение обратных задач методом Ньютона

Пусть необходимо найти минимум функции многих переменных f(X), где X = (x1, x2, x3, …, xn). Эта задача эквивалентна задаче нахождения значений X, при которых градиент функции f(X) равен нулю:

grad(f(X)) = 0.                                                    (11)

Применим к (11) метод Ньютона:

,        (12)

где j = 1, 2, 3, …, m – номер итерации; H(X) – гессиан функции f(X).

Напомним, что гессиан функции – это симметричная квадратичная форма, описывающая поведение функции во втором порядке:

где , функция f(X) задана на n-мерном пространстве вещественных чисел.

В более удобном для вычислений виде формулу (12) можно представить в виде

.

Метод Ньютона прост в реализации, но нахождение матрицы Гессе сопряжено с большими вычислительными затратами из-за необходимости нахождения большого количества частных производных.

Метод сопряженных градиентов

Метод сопряженных градиентов является методом нахождения локального минимума (максимума) функции с использованием информации о ее значениях и градиенте.

Метод сопряженных градиентов математи- чески существенно сложнее метода Ньютона. Не углубляясь в особенности метода, рассмотрим итерационный алгоритм нахождения экстремума функции: , где k – номер направления поиска; j – номер итерации, условие перпендикулярности направлений поиска k+1 и k: , направление поиска: , где .

Метод сопряженных градиентов математически сложнее метода Ньютона, но требует меньше вычислений.

Метод Ньютона и метод сопряженных элементов реализованы в надстройках Excel Microsoft.

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

Приведенная выше математическая модель, обобщенная на случай 3D, была описана в среде разработки Microsoft Visual Studio 15 на языке программирования C++. Для представления исходных данных, описывающих геолого-физические свойства пласта, был выбран формат в виде таблицы с использованием ключевых слов для обозначения блоков данных. Выбор представления обусловлен возможностью редактирования данных пользователем без использования специализированных программных средств (например в Microsoft Excel), а также экспорта данных из других симуляторов. При редактировании большого объема исходных данных пользователю необходимо хранить метаданные, относящиеся к ячейкам таблицы. Для реализации предложен алгоритм фильтрации, который удаляет метаданные, маркированные символами «//» из каждого блока при считывании файла. После успешного считывания исходных данных производятся проверки на целостность данных, размеры связанных блоков, размеры связанных интервалов перфораций и слоев, отсутствие пустых ячеек внутри таблицы. При успешном считывании и прохождении всех проверок пользователю выводится список выполненных операций. В случае ошибки пользователь также получает сообщение с указанием позиции в блоке данных и способах устранения ошибки.

На следующем этапе выполняется алгоритм моделирования заводнения нефтяных месторождений с применением многопоточности. Создаются отдельные потоки для расчета статических данных, а затем очередь обработки функций, рассчитывающих динамические данные. На данном этапе производится расчет изменяемого шага по времени; в целях сокращения используемой оперативной памяти хранятся данные только по трем шагам. Далее происходит формирование выходного файла, в котором присутствуют ключевые слова, описание размерности данных, диапазоны адресов ячеек данных для блоков, а также нумерация скважин и слоев. Предлагаемый формат выходных данных позволяет произвести их анализ и обработку (например в Microsoft Excel), в том числе построить диаграммы без дополнительной подготовки.

Созданная программа получила название Oil production calculator, или сокращенно OPC. Она была проверена с помощью сравнительного тестирования с симуляторами Dz10 Каца Р.М. и Волгина Е.Р., Rubis Kappa Engineering, Eclipse Schlumberger. Получено хорошее совмещение результатов расчетов с помощью ОРС и результатов, полученных с помощью вышеописанных программ.

Новый подход к решению задач оптимизации при моделировании разработки нефтяных месторождений

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

Рассмотрим простейшую секторную модель с 4 скважинами – 1 нагнетательной и 3 добывающими. Параметры модели следующие:

-     размеры модели в плане – 1 500´1 500 м;

-     абсолютная проницаемость – переменная (см. табл. 1);

-     пористость – переменная;

-     начальное пластовое давление – 339 атм.;

-     начальная водонасыщенность – 0,2 д.ед.;

-     объемный коэффициент воды – 1,02 м3/м3; вязкость воды – 0,36 сПз; сжимаемость воды – 4,7*10-5 1/атм.;

-     объемный коэффициент нефти – 1,55 м3/м3; вязкость нефти – 0,397 сПз; сжимаемость нефти – 1,01´10-4 1/атм.; растворимость газа в нефти – 203,5 м3/м3; давление насыщения нефти газом – 18,8 атм.;

-     сжимаемость породы – 4,7*10-5 1/атм.;

-     капиллярным давлением пренебрегается;

-     функции относительных фазовых проницаемостей (ОФП) приведены на рисунке 3;

-     дебит жидкости добывающих скважин – данные в таблице 1; минимальное забойное давление в добывающих скважинах – 18,8 атм.;

-     закачка воды в нагнетательную скважину без ограничений; забойное давление в нагнетательной скважине – данные в таблице 1;

-     скин-фактор для всех скважин – 0 ед.;

-     радиус скважин – 0,2 м;

-     срок разработки – 35 лет;

-     притока воды из-за контура нет.

Таблица 1

Параметры модели

Table 1

Model parameters

Примечание: 1 – нагнетательная скважина, 2, 3, 4 – добывающие.

Добывающие скважины управлялись дебитом, а нагнетательная – забойным давлением.

Расстановка скважин в плане и конфигурация суперэлементов показаны на рисунке 4.

Сведения о длине совместного ребра соседних ячеек для суперэлементной модели и расстоянии между центрами соседних ячеек (скважинами) приведены в таблице 2.

Таблица 2

Геометрические параметры суперэлементной модели

Table 2

Geometric parameters of a superelement model

Итак, в данном случае задача оптимизации состоит в том, чтобы уточнить значения проницаемости, используемые в модели, зная историю разработки. Имея оценку проницаемости по данным геофизических и гидродинамических исследований, а также исследований керна, с помощью генератора случайных чисел рассчитаем n вариантов проницаемости по 4 штуки в каждой группе k1, k2, k3, k4, где цифры 1–4 соответствуют номерам ячеек. Примем n = 20. Рассчитаем значения показателей модели для каждой группы значений проницаемости. Запишем целевую функцию следующим образом:

где Qo – накопленная добыча нефти по скважине; Qw – накопленная добыча воды по скважине; Qiw – накопленная закачка воды; индекс c – расчетные показатели; индекс f – фактические показатели; индекс i – номер ячейки. Для изучения эффективности предлагаемого метода в качестве фактических показателей примем показатели, полученные путем моделирования при значениях проницаемостей k1 = 70, k2 = 53, k3 = 44, k4 = 38 мД.

Таким образом, для каждого набора проницаемостей ячеек 1–4 получено значение целевой функции F. Из опыта анализа разработки известно, что накопленная добыча и закачка прямо пропорциональны проницаемости либо ее логарифму. Составим следующие уравнения для целевой функции F:

,                      (13)

               (14)

         (15)

где a1–a9 – постоянные коэффициенты.

С помощью регрессионного анализа найдем коэффициенты соотношений (13)–(15), используя полученную таблицу F от k1, k2, k3, k4. Затем, применив метод сопряженных градиентов к (13)–(15), определим проницаемости k1, k2, k3, k4, приняв начальное приближение проницаемости 50 мД (для такого маленького числа параметров (4 штуки) не важно, какой метод использовать – Ньютона или сопряженных градиентов). Сравним результаты определения проницаемости, полученные с помощью применения методов Ньютона и сопряженных градиентов непосредственно к модели фильтрации, а также полученные путем применения метода сопряженных градиентов к регрессионным моделям (13)–(15) (см. табл. 3). Начальное приближение проницаемости – 50 мД. Для оценки достоверности определения проницаемости введем ошибку по проницаемости следующим образом:

,                                      (16)

где индекс c – расчетная проницаемость; индекс f – фактическая проницаемость; индекс i – номер ячейки. Фактически соотношение (16) представляет среднюю относительную ошибку определения проницаемости.

Из таблицы 3 видно, что наиболее точные значения проницаемости можно получить, применяя метод сопряженных градиентов к модели фильтрации (ошибка – 11,8 %). Однако в случае, когда это невозможно (например, из-за большого числа иско- мых параметров), используя регрессионную модель и метод сопряженных градиентов (либо метод Ньютона), можно получить значения проницаемости с ошибкой не более 15,0 %. При этом важно правильно выбрать вид регрессионной модели.

Заключение

В работе описана суррогатная модель заводнения нефтяной залежи при пластовом давлении выше давления насыщения. Модель представляет собой полностью явную численную схему на базе концепции суперэлементов. Это позволяет реализовать данный подход для любого месторождения без использования специального ПО.

Рассмотрены подходы к решению задач оптимизации при суперэлементном моделировании разработки нефтяных месторождений. При этом описано использование метода Ньютона и метода сопряженных градиентов. В качестве альтернативы при большом количестве искомых параметров (и/или больших размерах модели) предложено применять методы теории оптимизации не непосредственно к модели фильтрации, а к регрессионной модели, построенной на базе экспериментов, проведенных на модели фильтрации. Рассмотрен пример уточнения значений проницаемости при адаптации модели с четырьмя скважинами. На этом примере показана приемлемая точность такого подхода. Следует подчеркнуть, что важно правильно выбрать вид регрессионной модели.

Работа выполнена при поддержке Программы фундаментальных исследований РАН № I.33П, НИР № 0065-2015-0111.

Литература

1.     Азиз Х., Сеттари Э. Математическое моделирование пластовых систем. М.–Ижевск: Изд-во Ин-та комп. исследований, 2004. 416 с.

2.     Каневская Р.Д. Математическое моделирование гидродинамических процессов разработки месторождений углеводородов. М.–Ижевск: Изд-во Ин-та комп. исследований, 2002. 140 с.

3.     Кац Р.М., Волгин Е.Р., Афанаскин И.В. Численное моделирование двухфазной фильтрации нефти и воды // Тр. НИИСИ РАН, 2014. Т. 4. Вып. 2. С. 141–148.

4.     Coats K.H., Thomas L.K., Pierson R.G. Compositional and black oil reservoir simulation. SPE 29111, 13th SPE Sympos. on Reservoir Simulation, 1995, pp. 149–162.

5.     Crichlow G.B. Modern reservoir engineering – a simulation approach. New Jersey: Prentice Hall. Inc., 1979, p. 354.

6.     Peaceman D.W. Fundamentals of numerical reservoir simulation. Amsterdam: Elsevier Scientific Publishing Co., 1977, p. 191.

7.     Хисамутдинов Н.И., Хасанов М.М., Телин А.Г. [и др.]. Разработка нефтяных месторождений. В 4 т. Т. 1. Разработка нефтяных месторождений на поздней стадии. М.: Изд-во ВНИИОЭНГ, 1994. 240 с.

8.     Мазо А.Б., Булыгин Д.В. Суперэлементы. Новый подход к моделированию разработки нефтяных месторождений // Нефть. Газ. Новации. 2011. Вып. 11. С. 6–8.

9.     Булыгин Д.В., Мазо А.Б., Поташев К.А., Калинин Е.И. Геолого-технические аспекты суперэлементной фильтрационной модели нефтяных месторождений // Георесурсы, 2013. Вып. 3. С. 31–35.

10.   Афанаскин И.В., Егоров А.А., Колеватов А.А. Экспресс-моделирование заводнения нефтяных месторождений с помощью концепции суперэлементов // Вестн. кибернетики. 2016. Вып. 2. С. 153–163.