Разработка и реализация численных методов моделирования многокомпонентной неизотермической фильтрации Шевченко Александр Валерьевич

Содержание к диссертации

Введение

1 Математические модели 10

1.1 Введение 10

1.2 Фазовый и компонентный состав 11

1.3 Уравнение баланса количества вещества 13

1.4 Уравнение баланса энергии 14

1.5 Замыкающие соотношения

1.5.1 Фазовое равновесие 14

1.5.2 Химические реакции 15

1.5.3 Свойства компонентов 17

1.5.4 Свойства фаз 18

1.5.5 Уравнения состояния фаз 18

1.5.6 Корреляции вязкости 21

1.5.7 Относительные фазовые проницаемости 21

1.5.8 Литотипы 22

1.5.9 Скважины 23

1.5.10 Начальные и граничные условия 24

1.6 Замкнутая система 24

2 Численный метод 25

2.1 Расчетные величины 25

2.2 Уравнение пьезопроводности 26

2.3 Расщепление по процессам

2.3.1 Перенос массы и энергии 29

2.3.2 Теплопроводность 33

2.3.3 Химические реакции 33

2.4 Блок фазового равновесия 34

2.4.1 Постановка задачи и допущения 37

2.4.2 Энтальпия 39

2.4.3 Метод решения 40

2.5 Аппроксимация пространственных величин з

2.6 Общая схема алгоритма 41

2.6.1 Инициализация 41

2.6.2 Шаг по времени 42

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

3.1 Верификация и исследование блока фазового равновесия 44

3.1.1 Верификация 44

3.1.2 Исследование алгоритма 45

3.1.3 Пример расчётной задачи

3.2 Верификация и исследование симулятора 51

3.3 Моделирование вибросейсмического воздействия

3.3.1 Постановка задачи 56

3.3.2 Результаты моделирования 59

Заключение 61

Список литературы 62

Список рисунков 65

• Уравнение баланса количества вещества
• Относительные фазовые проницаемости
• Химические реакции
• Верификация и исследование симулятора

Введение к работе

Актуальность и степень разработанности темы исследования

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

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

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

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

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

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

1. Изучить свойства уравнений неизотермической многокомпонентной фильтрации, определить функции, отвечающие за гиперболические и параболические свойства системы (если таковые имеются);

2. Разработать вычислительные алгоритмы для численного расчета фазового равновесия и кинетики химических реакций;

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

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

Основные принципы, положенные в основу построения вычислительного алгоритма, состояли в следующем:

4. Консервативность и монотонность разностной схемы для всех компонентов смеси и энергии;

5. Неявная аппроксимация всех переменных;

6. Расщепление по физическим процессам.

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

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

Цели и задачи работы Целями данной работы являются:

7. Разработка новых подходов к численному моделированию неизотермической многокомпонентной фильтрации, позволяющих учесть разложение/образование вещества скелета и динамическое изменение пористости и проницаемости вмещающей пористой среды;

8. Разработка вычислительных алгоритмов интегрирования уравнений, входящих в соответствующие гидродинамические модели фильтрации;

9. Реализация разработанных вычислительных алгоритмов в формате программного комплекса, его верификация и аппробация.

Задачи работы:

10. Разработать единый подход к численному моделированию многофазной фильтрации с динамическим учётом изменения пористости и проницаемости среды;

11. Конкретизировать выбранные модели, подобрав подходящие замыкающие соотношения для уравнений (корреляции вязкости, уравнения состояния фаз, зависимости относительных фазовых проницаемостей), условия на внешних границах расчётной области, а также модели скважин, учитывающие закачку и добычу многокомпонентной многофазной смеси;

12. Применить конечно-объёмную методологию для разработки устойчивого консервативного вычислительного алгоритма для решения основных уравнений модели неизотермической фильтрации (уравнения баланса количества вещества и энергии);

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

5. Реализовать разработанные алгоритмы в едином программном комплексе;

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

7. Провести моделирование тепловых МУН на керогеносодержащих месторождениях;

8. Провести моделирование ВСВ, реализованного на одном из опытных участков, и сравнить полученные результаты с историей разработки месторождения.

Научная новизна

Разработан подход, позволяющий описать процессы неизотермической фильтрации многокомпонентной смеси с фазовыми переходами и химическими реакциями. Разработан новый метод расчёта фазового равновесия многокомпонентной смеси, в котором единообразно учитываются подвижные и неподвижные фазы. Этот подход позволяет учитывать динамические изменения пористости и проницаемости коллектора. Предложенный метод основывается на задаче минимизации потенциала Гиббса системы «скелет-флюиды». Разработанные подходы и методы реализованы в виде вычислительного комплекса. Сравнение созданного вычислительного комплекса с существующими аналогами показало согласование результатов расчётов на спектре общих задач. Впервые для моделирования пилотных разработок опытных участков на реальных месторождениях был применён программный комплекс, способный учитывать динамические изменения пористости и проницаемости коллектора. Результаты моделирования воспроизвели историю эксплуатации с достаточной точностью.

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

изменения пористости и проницаемости скелета. На основе этих подходов создан программный комплекс для моделирования указанного класса задач.

Сравнение результатов расчётов с аналитическими решениями и результатами коммерческих симуляторов подтверждает корректность предложенных подходов и разработанных численных алгоритмов. Созданный программный комплекс можно применять для расчёта эффективности проектов разработки месторождений трудноизвлекаемых углеводородов.

Разработанная модель пиролиза керогена и подход динамического учёта изменений пористости и проницаемости коллектора были положены в основу расчётов по исследованию эффективности тепловых МУН на месторождениях высоковязких нефтей. Выполнена реконструкция истории разработки опытного участка Жирновского месторождения, на котором осуществлялось вибросейсмическое воздействие. Использование прямого численного моделирования на разработанном программном комплексе способно существенно повысить достоверность ожидаемых результатов при различных схемах разработки месторождения и, как следствие, позволит выбрать оптимальную технологию.

Работа поддержана рядом государственных грантов:

Грант ФЦП «Научные и научно-педагогические кадры инновационной России» «Механизмы вибросейсмического воздействия, повышающие нефтеотдачу пласта и интенсифицирующие добычу нефти» 2010-2012 гг.

Разработка программного комплекса является одним из проектов в лаборатории флюидодинамики и сейсмоакустики МФТИ (финансирование из средств Программы повышения конкурентноспособности МФТИ «5топ100»)

Грант РНФ №15-11-00015 «Математическое моделирование физико-химических воздействий на пласт и скважинных процессов для целей оперативного управления добычей тяжёлой нефти»

Методология и методы исследования

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

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

Положения, выносимые на защиту

Положения, выносимые на защиту, отражены в основных результатах работы, приведенных в конце автореферата.

Степень достоверности и аппробации работы.

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

Результаты диссертации опубликованы в 13 работах, из которых 3 в изданиях, рекомендованных ВАК РФ [–], а одна содержится в международной базе SPE OnePetro [].

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

1. 53 — 57 научной конференции МФТИ, МФТИ, Долгопрудный, 2010-2014 (6 докладов);

2. 3-rd international conference «High Performance Computing» (HPC-UA’13, HPC Day’13), КПИ, Киев, 2013;

3. SPE Russian Technology Conference, Moscow, 2015.

В рамках работы было получено свидетельство на программу для ЭВМ [].

Уравнение баланса количества вещества

В уравнениях (1.1), (1.2), а также уравнениях состояния фаз, о которых речь пойдёт ниже, используются массовая и молярная плотности компонентов. В общем случае, их зависимости от давления р и температуры Т могут быть достаточно сложными (см. раздел 1.5.4). Для корректного пересчёта одной величины в другую необходимо знать молярную массу Mj компонента. Эта величина является постоянной и не зависит от внешних условий.

В силу аддитивности, молярная энтальпия фазы а, фигурирующая в уравнении (1.5) имеет вид: где hia = Ыа(р,Т) — молярная энтальпия компонента г в фазе а, в общем случае зависящая от давления р и температуры Т.

В соответствии с этим, для компонентов, которые содержатся только в одной фазе (тяжелая нефть, нерастворимые газы, неподвижные компоненты) требуется задать только одно соотношение для энтальпии; в то время как для испаримых (лёгкие нефти и вода) — два. В симуляторе принята следующая форма зависимости молярной энтальпии от давления и температуры [Принятая также и Eclipse]: h = А(Т — Tref) -\ (Т — Tref) , Ср = А + В(Т — Tref) При этом параметры для энтальпии неподвижных компонент получаются из объёмной теплоёмкости С следующим образом: Р Также, для испаримых компонентов задаётся теплота фазового перехода жидкость-газ (для компонента НгО — вода газ) [Принятая в CMG STARS + аналогичная формула в Eclipse]:

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

Уравнение состояния (УС) фазы необходимо для расчёта объёма, который она занимает, при заданном давлении и температуре [Ссылка!]. В симуляторе возможно задание четырех видов уравнений состояния фаз: “Compressible” — УС, учитывающее сжимаемость и термическое расширение [Ссылка!]; “Ideal” — УС, описывающее идеальный газ с фиксированным Z-фактором [Ссылка!]; “WaterEos” — УС, описывающее поведение воды в пластовых условиях; “PengRobinson” — УС Пенга-Робинсона [Ссылка!].

Иногда вместо использования УС фазы, считают, что компоненты в фазе составляют идеальный раствор (объём фазы является суммой объёмов компонентов) [Ссылка! Документация STARS]. В таком случае требуется задание УС всех компонентов, входящих в состав фазы.

Уравнение состояния воды Данное УС может использоваться только для НгО и является специфичным для разработанного симулятора.

Такое УС позволяет учесть сжимаемость Н20 в пластовых условиях и считать этот компонент несжимаемым на поверхности. Использование УС на поверхности обусловленно пересчётом дебитов добывающих скважин из молярных величин в объёмные.

Уравнение состояния с постоянными коэффициентами сжимаемости и теплового расширения

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

В этом уравнении состояния КІ(ШІ) имеет разрыв порядка 0.1% в точке 0.49. Возможно, это вызвано недостаточной точностью числовых констант, либо опечатками в основных работах по этой теме [Ссылка!]. Случай больших ацентрических факторов маловероятен для рассматриваемых задач и данная погрешность не повлияла ни на один расчёт в рамках этой работы.

В уравнение (1.2) входит ещё один фазовый параметр — относительная проницаемость фазы [Ссылка!]. Существует несколько вариантов задания этой функции: табличное задание, степенной закон, “Eclipse” (вариант, который используется по умолчанию в симуляторе Eclipse) и корреляция согласно первой модели Стоуна. Стоит отметить, что часто именно способ вычисления относительной фазовой проницаемости может существенным образом повлиять на процесс фильтрации и, как следствие, требуется тщательное изучение кернов месторождений и их фильтрационных свойств, для получения достоверных результатов при моделировании. Подробный обзор различных моделей относительной фазовой проницаемости можно найти в [21]. Степенной закон

Данная модель для вычисления относительной проницаемости нефти была предложена Стоуном в [22]. Основной особенностью данной модели является учёт эффекта повышения относительной проницаемости нефти в случае присутствия небольшого количества газа. В силу громоздкости формул, которые используются в данной модели, её описание приведено в приложении B. Более подробно о моделях относительных фазовых проницаемостей можно прочитать в [21].

Под литотипами (литологическими типами) обычно понимают типы породы, идентифицируемые по набору литологических признаков [Ссылка!]. В разработанном симуля-торе предполагается, что литотип задаётся пористостью , насыщенностями неподвижных фаз Si (которых формально может быть больше, чем одна), а также тензором абсолютной проницаемости К. При этом считается, что К является диагональным, то есть величины KXZ = 1KZX = 0. В случае, если в литотипе присутствует неподвижный компонент, который может входить в реакции с разложением на неподвижные и подвижные фазы, требуется задание дополнительной проницаемости, которая возникает при разложении этого компонента.

В симуляторе рассматривается несколько видов скважин:

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

Добывающие. Моделируются как источниковые члены в уравнениях баланса количества вещества и энергии. На перфорированных участках добывающих скважин задается забойное давление. Дебит по фазам, дебит компонентов и ток энергии вычисляется с учетом состояния смеси в блоке, содержащем соответствующий перфорированный участок. На добывающих скважинах также рассчитывается объёмный и массовый дебит компонентов на поверхности (при нормальных условиях).

Нагреватели Моделируются как источниковые члены в уравнении баланса энергии. При этом задаётся предельная мощность и максимальная температура, при которой работает нагреватель.

Математическая модель нагнетательных и добывающих скважин основана на формуле Писмана [Ссылка!], которая описывает приток фазы в прямоугольный блок. Эта формула выражает закон Дарси, но с безразмерными множителями, учитывающими распределение давления в окрестности скважины. Пусть ось z направлена вдоль скважины, тогда объем фазы а, поступающий в блок каждую секунду равен:

Относительные фазовые проницаемости

В данной работе, однако, считается, что часть компонентов не содержится в некоторых фазах (см. раздел 1.2). В этом смысле рассматривается ограниченная задача о фазовом равновесии. Компоненты смеси разделяются на два типа: активные — те, которые участвуют в установлении фазового равновесия, и инертные, роль которых ограничивается участием в балансе энергии. Принимается что активные компоненты могут образовать три фазы: жидкую (нефтяную), газовую и водную. Жидкая фаза содержит компоненты, не присутствующие в других фазах (тяжелые нефти) и компоненты, присутствующие в газовой фазе (легкие нефти). Газовая фаза наряду с легкими нефтями содержит нерастворимые газы (содержащиеся только в газовой фазе) и пары Н20. Водная фаза состоит только из одного компонента Н20. Количество фаз, на которое расслаивается смесь, заранее неизвестно и определяется термобарическими условиями и составом смеси. Инертные компоненты образуют скелет — набор неподвижных фаз.

Задача о фазовом равновесии, которая рассматривается в данной работе, состоит в следующем. Заданы относительные молярные концентрации всех компонентов смеси Q, её молярная энтальпия г/ и давление р. Требуется определить на какие фазы расслоится смесь, их молярные доли и составы, температуру, а также производные dV/dci,dV/dr].

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

Задача минимизации потенциала Гиббса (при фиксированных температуре, давлении и составе смеси) ставится в пространстве независимых степеней свободы: молярных долей (относительно полного количества молей) Н20 в водной фазе и легких нефтей в жидкой фазе. Указанные величины однозначно описывают состояние системы (молярные доли и компонентный состав фаз). Условия минимизации потенциала Гиббса в такой постановке дают критерии фазового равновесия для сосуществующих фаз, а балансы количества молей компонентов (уже учтенные в аналитическом виде потенциала) приводят к уравнениям Речфорда–Райса для молярных долей фаз. Множество, на котором происходит минимизация ограничено и описывается набором естественных условий неотрицательности молярных долей компонентов.

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

При вычислении производных возникает ряд дополнительных проблем. Во-первых, производные перестают существовать на границе области минимизации. Во-вторых, при вариации параметров (температуры, давления, состава) может измениться структура множества ограничений. Эта проблема возникает в случае, когда какие-то компоненты смеси отсутствуют.

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

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

В результате расчета фазового равновесия в п.5 предлагаемого алгоритма активные компоненты смеси всегда образуют трехфазную систему. Если система фактически (т.е. с исходным модельным потенциалом Гиббса) находится в одно- или двухфазном состоянии, то молярные доли "паразитных"фаз не должны превышать некоторого порога, регулируемого высотой внутреннего барьера.

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

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

Рассматривается многокомпонентная система, которая может расслаиваться на несколько фаз. Часть из компонентов (инертные) не участвует в фазовых превращениях и формирует инертные фазы. Другие компоненты (активные) могут потенциально образовывать три подвижные фазы: жидкую (нефтяную) фазу (индекс ), газовую фазу (индекс ) и водную фазу (индекс ).

Химические реакции

Основной задачей моделирования было качественное повторение эффекта увеличения нефтеотдачи при проведении ВСВ. По данным полевых и лабораторных экспериментов, при проведении ВСВ отмечаются факты дегазации (выделения газа из нефти) с объёмными долями газа до 8% при продолжительном ВСВ. Предположительно, такое поведение связано с тем, что давление в резервуаре поддерживалось на несколько атмосфер выше давления раз-газирования, и возмущения в резервуаре, индуцируемые сейсмическими волнами, позволяют выделиться газу. При этом обратный процесс растворения газа протекает медленно, что позволяет рассматривать газ как отдельную фильтрующуюся фазу на временах сравнимых с опытным периодом (порядка одного года).

В рамках лабораторных исследованих были построены зависимости относительных про-ницаемостей для трехфазной системы Нефть-Вода-Газ (см. рис. 3.16). Приведенные данные свидетельствуют о сложной зависимости относительной проницаемости нефтяной фазы от насыщенностей всех фаз (а не только от насыщенности нефтяной фазы). Можно видеть, что для месторождения с высокой обводнённостью (левый нижний угол диаграммы), возможны ситуации, при которых нефть, обладающая нулевой проницаемостью, станет подвижной, при добавлении газа в флюид, что приведёт к дополнительной добыче нефти. Подходящей корреляцией относительных фазовых проницаемостей для данного случая является модель Стоуна (см. приложение B).

Помимо особенного вида фазовых проницаемостей, большое значение имеет повышение давления в резервуаре за счёт выделения газа (ранее растворённого в нефти). Повышение давления влиет на скорости фильтрации, увеличивая общее количество добываемой жидкости.

Во время проведения ВСВ с помощью скважинных приборов были измерены величины акустического шума в пласте. Акустический шум может являться индикатором энергонасыщения породы, и, как следствие для данного конкретного месторождения, количества газа, вышедшего из нефти. Так, на рис. 3.17 красным пунктиром отмечена область, в которой наблюдался рост акустического шума. Зелёная линия соответствует выбранному расчётному профилю. Данный профиль приближен к срезу месторождения по которому в отчёте были приведены данные о конфигурации пласта. На всём протяжении этой линии пласт остаётся практически горизонтальным, за исключением области вблизи нагнетательной скважины 268, где пласт начинает уходить вглубь (этим фактов пренебрегаем при численном моделировании). В расчёте рассматриваются три добывающие скважины, находящиеся в непосредственной близости от расчётного профиля — это скважины 915, 914 и 182; а также две закачивающие — 268 и фиктивная скважина, образованная суммарным действием скважин 153 и 60. На нагнетательных скважинах задается расход воды в соответствии с установленным режимом работы на месторождении (расход на фиктивной скважине равен сумме расходов на скважинах 153 и 60). Стоит отметить, что в силу использования 2D симулятора, прямого соответствия результатов моделирования с полевыми данными получить нельзя, поскольку на добычу из рассматриваемых скважин влияют также соседние добывающие (177, 913, 101) и нагнетательные скважины (125).

На рисунке 3.18 изображено распределение абсолютных проницаемостей, используемое в расчёте. Более светлые области соответствуют большим значениям проницаемости. Области чёрного цвета обозначают литотипы с нулевой проницаемостью. Данное распределение было восстановлено из данных об абсолютной проницаемости, полученных в добывающих скважинах. Видна сложная структура коллектора, в данном случае разделённого на два независимых горизонта. В отчёте присутствуют также данные профиля притока флюида вдоль добывающий скважин. Для удовлетворения этих зависимостей были добавлены три области с нулевой проницаемостью, через которые проходят скважины. Для скважин 182 и 914 эти области — естественные образования, в то время как основание скважины 915 было зацементировано.

В отчёте также присутствуют данные о текущей (на момент начала проведения ВСВ) водо- и нефтенасыщенности пластов в добывающих скважинах. Эти данные стали не только начальными данными для расчёта, но из них также были вычислены значения остаточных

Распределение абсолютных проницаемостей вдоль линии среза. Более светлые области соответствуют большим значениям проницаемости. нефте- и водонасыщенностей. В силу того, что добыча газа на забое наблюдалась только на одной скважине, остаточная насыщенность газа была положена 0.1 — значение большее мак 59 симальной насыщенности газа, которая может выделиться при проведении ВСВ в условиях месторождения (0.08).

Область вибрационного воздействия на рисунке выделена красной рамкой. Виброисточник моделировался как дополнительная скважина, которая добавляет в область на каждом временном шаге моль/м3 газа, где — так называемый генерационный потенциал нефти — количество газа, которое выделяется из 1 моля нефти за единицу времени. Генерационный потенциал подбирается таким образом, чтобы за время работы вибросейсмического источника (60 дней) выделилось количество газа, имеющего объёмную долю примерно равную 0.08 для ячейки, полностью насыщенной нефтью.

На рисунке 3.19 изображено сравнение накопленной добычи нефти в скважине 915 без (красный график) и с (зелёный график) ВСВ. Наблюдается три характерных периода при разработке: период проведения ВСВ, область увеличения добычи и период выхода на первоначальный режим разработки. Увеличение добычи на месторождении наблюдалось в течение одного года, что соответствует данным расчёта.

Верификация и исследование симулятора

В точке 2.865 наблюдаются нефизические осцилляции погрешности при достаточно малых є. Это явление вызвано численной ошибкой, которая накапливается при вычислениях с большой точностью, поскольку указанная точка находится в окрестности критической для чистого вещества С5Н12. Параметры критической точки используются при вычислении корреляции энтальпии парообразования (разности энтальпий между фазами) [27]. Для того чтобы избежать разрыва производной в данной корреляции, она также сглаживается с погрешностью порядка 10-6. Тем не менее, нефизические осцилляции не наблюдаются на других графиках (зависимостях долей СбН12 в жидкой и Н20 в водной фазах от энтальпии) и на данном графике вплоть до точности є = 10-6. Всё это свидетельствует о том, что на их появление оказала влияние комбинация факторов, указанных выше.

Отметим, что ошибка решения ведёт себя линейно в зависимости от погрешности практически во всей расчётной области, однако в точках изменения количества фаз точность падает вплоть до у/є. Таким образом, для обеспечения одинаковой погрешности во всём диапазоне температур, стоит вести расчёты в областях смены количества фаз с повышенной точностью.

C10H22 в жидкой фазе, зелёный —- доле C5H12 в жидкой фазе и синий –– доле H2O в водной фазе. Можно видеть, что точки-узлы наблюдаются и для других переменных: для C5H12 имеются две точки-узла в трёхфазной области и одна —- в двухфазной (водная фаза отсутствует); а для H2O одна точка в трёхфазной области. Также отметим, что на графиках, соответствующим C5H12 и H2O , отсутствуют осцилляции в окрестности точки 2.865.

На рисунке 3.6 представлен сводный график работы алгоритма на модельной задаче. Синяя, красная и зелёная линии представляют молярные доли компонентов в жидкой и водной фазах. Розовая линия является афинно преобразованным графиком зависимости температуры от энтальпии, а чёрная линия представляет собой количество итераций метода Ньютона, затраченное для нахождения точки равновесия. Последнюю зависимость проанализируем отдельно. Как и следовало ожидать, в областях смены количества фаз увеличивается количество итераций, которое необходимо для нахождения равновесия. При этом количество итераций возрастает незначительно —- максимум на 10-12. В целом, алгоритм быстрее схо 49

дится в трёхфазной области (0.288 – 1.934) и медленнее в областях, где присутствует газовая фаза. В области низких энергий (отсутствует газовая фаза) в некоторых точках при больших точностях ( 10-9) наблюдается отсутствие сходимости и метод прекращает работу из-за достижения максимально допустимого количество итераций (в модельной задаче использовалось значение 100, на графике изображенное величиной 0.5). В этой области метод достаточно быстро достигает окрестности точки равновесия, после чего долгое время итерации происходят в ней. Отсутствие сходимости вызвано накоплением численных ошибок, вызванных большой обусловленностью матриц, возникающих в методе Ньютона. В реальных задачах такая точность расчётов является чрезмерной (достаточная точность = 10-6). При использовании меньшей точности подобные проблемы возникать не будут. Время исполнения одной итерации метода было примерно одинаковым во всём диапазоне температур и составляет несколько микросекунд (в зависимости от оборудования, на котором производятся расчёты).

Разработанный алгоритм был протестирован при расчёте фазового равновесия в пластовых условиях месторождения А. Данное месторождение характеризуется наличием нескольких инертных компонентов, составляющих скелет, двумя газами CO2 и N2 в незначительных количествах, растворённых в нефти в пластовых условиях; а также несколькими фракциями нефти, объединёнными в четыре группы, условно названные C1, C2, C10 и C20. Считается, что первые три фракции могут испаряться, а последняя является тяжёлой нефтью. Режим разработки месторождения предполагает тепловое воздействие, что приводит к необходимости расчёта фазового равновесия, как составляющего компонента гидродинамического симуля-тора. Пластовое давление составляет порядка 100 атмосфер, в то время как температура изменяется от 100 до 300 градусов Цельсия. Все расчёты проводились с точностью = 10-6, растворимые газы моделировались как лёгкие нефти, корреляции для всех зависимостей (энтальпий фаз, констант фазового равновесия и энтальпий парообразования) взяты из открытых источников [27].

На рисунке 3.7 показана зависимость долей компонентов в жидкой (для лёгких нефтей) и водной (для H2O) фазах от температуры (в кельвинах). Цветами обозначены компоненты: синий H2O, красный CO2, зелёный N2, фиолетовый C1, оранжевый C2 и коричневый C10. Начиная с температуры порядка 282 К, наблюдается дегазация нефти и испарение её лёгких фракций и H2O. Полное испарение происходит в рассматриваемой системе вблизи 490 K, т.е. 225 C. Из графика можно заключить, что испарением фракции C10 в целом можно пренебречь, поскольку оно становится значимым только в областях с температурами большими 500 K. вблизи точек смены числа фаз, на которое расслоится смесь. Среднее количество итераций, необходимое для достижения заданной точности в трёхфазной области составляет 7-10, а в двухфазной —- 13-15. При этом максимальное количество итераций, требуемое только в узкой окрестности точек смены количества фаз —- 23. За исключением этих областей алгоритм устойчиво получает решение за 18 и менее итераций. Вышесказанное свидетельствует о его надёжности и возможности равномерно распределять нагрузку (что не так для «стандартного» алгоритма расчёта фазового равновесия). 3.2 Верификация и исследование симулятора

Для обеспечения адекватности моделирования была проведена предварительная проверка разработанного симулятора и его сравнение с коммерческим симулятором Schlumberger Eclipse 300.

В верификационном расчётном сценарии рассматривался трёхслойный коллектор, центральный слой которого был более проницаемым. Нефтяная фаза состояла из четырёх компонентов, три их которых являлись лёгкими нефтями, а один — тяжелой нефтью. При этом, в пластовых условиях в нефти были также растворены два газа. PVT-свойства флюида, закачиваемой воды и газов были взяты из открытых источников. В области имеются две скважины: нагнетающая и откачивающая. Обе описывают приток вещества в соответствии с формулой Писмана. Добывающая скважина осуществляет работу при фиксированном давлении на устье, в то время как на закачивающей указаны как максимальное давление, так и максимальный поток закачиваемой воды. Расчёт вёлся с прогнозом на 10 лет.

На рисунке 3.9 приведено сравнение полей насыщенности воды в определённый момент времени. Наблюдается схожая картина фильтрационного течения: центральный слой был практически полностью вытеснен (до значения критической нефтенасыщенности), два менее проницаемых пластах только частично вытеснены. Не смотря на одинаковые характеристики менее проницаемых пластов, фильтрационная картина в них разная за счёт учёта гравитационных сил. Отметим, что в результатах, полученных с помощью Schlumberger Eclipse наблюдается немонотонное поведение решений в областях, где насыщенности близки к критическим значениям.

На рисунке 3.10 приведено сравнение дебитов добытой нефти. Наблюдается практически полное согласование графиков накопленной добычи до момента прорыва центрального слоя, после чего разработанный симулятор BSD демонстрирует меньший темп добычи по сравнению с Eclipse. Разница в накопленной добыче не превосходит 4%. Аналогичные результаты демонстрируются и на графике добычи газа.