В условиях уменьшения запасов легких нефтей возникает потребность в разработке месторождений высоковязких нефтей и природных битумов. Однако для разработки таких ресурсов нужны специальные технологии. Как правило, высоковязкие нефти разрабатываются с использованием тепловых методов. Одним из перспективных тепловых методов является высокочастотное электромагнитное воздействие [1–6]. Одним из сдерживающих факторов при внедрении данного метода является отсутствие инструментов для моделирования воздействия.
Высокочастотное электромагнитное воздействие на пласт и призабойную зону основано на использовании особенностей термогидродинамических процессов в продуктивных коллекторах [7–9]. Поле действует, в частности, на компоненты тяжелых нефтей, образующие сложные конгломераты и сетчатые структуры с полярными элементами, такие как смолы и асфальтены, потому что они являются полярными соединениями. Под действием внешнего поля полярные составляющие соединений начинают ориентироваться вдоль направления поля, а так как это поле меняет свое направление с высокой частотой, то происходит разрушение этих структур и выделение теплоты. Усиление воздействия ВЧ ЭМ нагрева проводится закачкой смешивающимся с пластовым флюидом агента, которое происходит за счет дополнительного уменьшения вязкости углеводородного флюида в результате смешивания и вследствие уменьшения требуемой температуры нагрева флюида [10–12].
Классические подходы к решению задачи фильтрации, в частности, при моделировании ВЧ ЭМ воздействия [13] используют модели простой сжимаемости, в решениях используются осредненные свойства флюидов и параметры гипотетических пластов, поэтому они не учитывают эффекты ограниченного растворения компонентов друг в друге, а при осуществлении расчета многофазных течений с фазовыми переходами такой подход некорректен. Для моделирования ВЧ ЭМ прогрева пласта на данный момент не существует инструмента, сочетающего возможности учета температурных изменений в пласте со смесимостью фаз и вычислительной простотой метода.
В данной работе осуществлена связь классической математической модели воздействия ВЧ ЭМ поля с реальными пластовыми условиями. Основная сложность задачи заключается в исключении описания PVT-свойств флюидов через сложные уравнения состояния многокомпонентной многофазной системы, которой является нефтяной пласт, с учетом того, что термобарические условия в процессе ВЧ воздействия и фильтрации флюидов будут непрерывно меняться.
Чтобы учесть то, что происходит в пласте в реальности, по нашему мнению, самый лучший путь — воспользоваться свойствами конкретного пласта напрямую, т.е. использовать результаты экспериментальных исследований пластовых нефтей месторождения, с помощью специализированных симуляторов создать модель свойств (PVT-свойств), согласованных между собой, настроить модель на эти результаты исследований, таким образом получить распределение свойств по реальному пласту в виде набора свойств флюидов. И эту базу данных (БД) свойств — плотности, объемных коэффициентов в зависимости от давления, температуры и насыщенностей (т.н. PVT-таблицы) в данной точке пласта — уже можно использовать в аналитической задаче для наблюдения поведения пласта в условиях воздействия ВЧ ЭМП, т.е. условиях изменения температуры, давления, концентраций компонент в потоке. БД будет использоваться вместо уравнений состояния системы, которые задаются в аналитических задачах.
Постановка данной задачи и краткое описание методики расчета свойств фаз, участвующих в совместной фильтрации при воздействии ВЧ ЭМ поля, приведена в работе [14]. В данной работе описаны подходы к решению задачи, инструменты для решения и тестирование задачи, приведены результаты расчетов оценки влияния внедрения нового метода расчета свойств флюида, произведены сравнительные расчеты двух моделей задачи: с использованием осредненных свойств флюидов и параметрами гипотетических пластов и с использованием реальных свойств и параметров.

В рамках описания физической постановки задачи рассмотрим процесс фильтрации в пористом пространстве многокомпонентного многофазного флюида, в частности, высоковязкой нефти с включениями, асфальтосмолистых и парафиновых отложений, смешивающегося растворителя, нагнетаемого в пласт. При воздействии высокочастотным электромагнитным излучением за счет поглощения энергии поля максимальный нагрев пласта и флюидов происходит не на стенке скважины, а на некотором расстоянии от нее в пласте, а также закачиваемый агент несет нагретую жидкость из призабойной зоны скважины в глубь пласта. Температура в рассматриваемой системе постепенно увеличивается и в некоторой области достигает температуры плавления твердых включений, где начинается фазовый переход твердых частиц в жидкое состояние. Расплав твердой фазы также вовлекается в процесс фильтрации закачиваемого растворителя и нефти. Вследствие дальнейшего нагрева среды ВЧ ЭМ полем область фазового перехода будет распространяться и в глубь пласта, и к забою скважины.
Для математического описания этого неизотермического процесса, по сути многофазной многокомпонентной неизотермической фильтрации флюидов в пористой среде, справедливы следующие фундаментальные законы и уравнения: закон сохранения вещества для каждого компонента (1), закон сохранения энергии (2), модель распределенных источников тепла [15, 16] при ВЧ ЭМ воздействии на диэлектрическую среду (3).

Начальные и граничные условия замыкают эту систему уравнений:

Здесь m — пористость, д.е., N_c — мольная концентрация, моль/м³; Q_s — количество вещества, привнесенного извне, моль/с; x^p_c — доля компонента с в фазе p; ξ^d_p — мольная плотность фазы, моль/м³; u→ поток флюида, м/с; v^→ — скорость изменения концентрации, моль/м²/с; A — площадь поверхности объема, м²; V — объем, м³; t — время, сек; Σf = (-x^p_c ξ^d_pu)— суммарный поток компонента, распределенный по выделенным фазам, моль/м²/с, C_c — эффективная теплоемкость ячейки, Дж/м³/K; C_r — мольная теплоемкость горной породы, Дж/моль/K; λ^t_c — эффективная теплопроводность, Вт/м²/К; ρ_r — мольная плотность горной породы, моль/м³; k_f — коэффициент, характеризующий долю площади занимаемого флюидом в перпендикулярном направлении передачи тепла, д.е.; λ^f_hf — теплопроводность флюида, Вт/м²/К; λ^r_hf — теплопроводность горной породы, Вт/м²/К; u^e_f — количество вещества, переносящего энергию, Вт/м²/К; C_p — теплоемкость фазы, Дж/моль/К, Q_f — температурный эффект фазовых переходов, Дж/с; Т_пл — температура плавления, K; ∆N_f — количество вещества, совершившего фазовый переход, моль/с; l_f — мольная удельная энергия фазового перехода, Дж/моль; H^(m)_v — функция Ханкеля; tg(δ) — тангенс угла диэлектрических потерь; N₀ — мощность излучателя, Вт; ω — круговая частота, Гц; ɛ и μ — относительная комплексная диэлектрическая и магнитная проницаемости среды; h — толщина пласта, м; r_wb — радиус скважины, м.

Согласно выражению (1), изменение количества вещества происходит за счет движения флюидов, за счет диффузии и притока или оттока, связанных с работой скважины. Закон сохранения энергии (2) показывает, что изменение энергии происходит за счет теплопроводности, конвективного переноса, воздействия ВЧ ЭМ поля, фазовых переходов. Энергия электромагнитного поля (3) частично поглощается диэлектрической средой, в результате чего в среде появляются распределенные источники тепла с объемной плотностью теплового потока.
Для численного решения системы уравнений (1–4) необходимо перейти к дискретной постановке задачи, после чего произвести линеаризацию. Ниже приведена система уравнений (1–4) в результате дискретизации и линеаризации (5):

Полученная система уравнений при записи для каждой ячейки решается численно с использованием метода Ньютона. Схема, по которой осуществлялось решение, приведена на рисунке 2.

На данном этапе основной проблемой является корректный учет плотности фазы, доли компонента в фазе, проницаемости и вязкости при расчете уравнения сохранения вещества для каждого компонента, т.е. моделирование свойств флюидов в пласте.
Основной проблемой при численном решении системы уравнений является способ расчета коэффициента проводимости λ_(c/b) (P,T,N) [17]

Для его нахождения необходимо рассчитать два очень важных параметра: плотность и вязкость, которые в отличие от классических задач являются функциями давления, температуры и концентрации компонентов.
На начальном этапе развития подземной гидромеханики при расчете проводимости пренебрегали влиянием сжимаемости [18] и изменением вязкости от давления, что позволяло решать подобную задачу аналитически. Дальнейшее развитие подземной гидромеханики было связано с численным решением сложнейших дифференциальных уравнений, это позволило учитывать более тонкие эффекты. На данном этапе стало возможным использовать линейные зависимости плотности (объемные коэффициенты) и вязкости от давления и при доработке подхода от температуры. Такой подход приемлем в случае расчета пластов, насыщенных только нефтью с низким давлением насыщения. Для пластов, насыщенных нефтью в контакте с газом, или в случае близости пластового давления и давления насыщения необходимо адаптировать подход расчета. В этих целях создана схема (PVT-модель), которую в специализированных симуляторах определяют через ключевые слова PVTO/PVTG. В рамках данной схемы можно выделить два элемента: ветвь насыщения, когда любое изменение давления сопровождается либо поглощением второй фазы, либо ее высвобождением, и вторая часть, в которую входят множество ветвей недонасыщенной нефти. Как правило, такой подход используют, когда вторая фаза — это газ, однако он принципиально применим и для смешения с растворителем. Преимуществами данного подхода для вычислений являются относительная простота и абстрагированность от физического наполнения.
Параллельно с усложнением подходов в подземной гидромеханике развивалась и термодинамика. Были предложены инструменты для построения моделей черной/летучей нефти гидродинамических симуляторах и физически наполненный подход моделирования PVT-соотношений с использованием сложных уравнений состояния (в рамках моделирования пластовых систем — композиционная модель флюида). У такого подхода есть существенный плюс: после построения данной модели и проверки ее на адекватность в интересующих пределах можно не задумываться о физическом содержании фазовых переходов при различных давлениях и температурах, что особенно важно для высокочастотного электромагнитного воздействия. Однако этот метод имеет и минусы: данный подход требует решения сложной системы уравнений, что приводит к значительному увеличению времени расчета модели, и для их построения требуются значительные трудозатраты высококвалифицированных специалистов.
Для того чтобы построить модель летучей нефти, можно взять конкретную пробу, детально исследовать процессы разгазирования и сжимаемости на ступенях разгазирования, или же исследовать множество различных проб только в области разгазирования, и на основе этих данных в специализированных PVT-симуляторах построить общую композиционную модель, на основе которой построить модель летучей нефти.
Если рассматривать абстрактно переход от одной модели к другой, то необходимо указать на то, что композиционная модель на основе множества коэффициентов, описывающих каждый компонент специального кубического уравнения, определяет некоторое термодинамическое фазовое пространство. А модель летучей нефти представляет собой некоторый срез этого пространства в несколько измененных координатах. Можно построить такой же срез для различных температур в диапазоне, который характерен для высокочастотного воздействия. Если графически обратно сложить эти срезы, мы получим фазовое пространство меньшей детальности, только теперь оно описывается не кубическим уравнением, а набором PVTO/PVTG подобных таблиц в измененных координатах. Таким образом, получен способ записи данных, достаточных для описания поведения нефти в контакте с растворителем под действием высокочастотного электромагнитного поля. Назовем данную модель термолетучей.
Отметим, что при переходе от композиционной модели к термолетучей разница между компонентом и фазой практически стирается. Это обусловлено тем, что каждый компонент может образовывать свою фазу, состоящую из себя самой и второго компонента, который также образует фазу.

Для использования таблиц с PVT-свойствами необходим также способ расчета свойств двух фаз, находящихся в контакте при данном давлении и температуре. Можно разбить нашу модель на две подмодели: состояние насыщенного флюида и ненасыщенного флюида. Первая, в отличие от изотермического случая, будет иметь вид не кривой, а плоскости, и для ее описания можно использовать двумерную линейную интерполяцию, где переменными будут давление и температура, а определяемыми параметрами — плотность, вязкость, растворимость одного компонента в другом. Эта подмодель описывает ситуацию, когда при заданных давлении и температуре заданный объем второй фазы не растворяется в первой, и они находятся в контакте. Вторая подмодель описывает фазу, когда второго компонента недостаточно для образования второй фазы. Для описания плотности и вязкости необходимо использовать метод n-мерной линейной интерполяции, где переменными будут давление, температура и содержание растворенного компонента.
Для определения свойств двухкомпонентной системы с заданной массой каждого компонента необходимо вначале определить на основе насыщенной подмодели предельное содержание второго компонента в первой фазе. Если максимальное содержание компонента в первой фазе меньше заданной массы, то для расчета плотности и вязкости первой фазы используется первая подмодель. Свойства второй определяются по аналогии с первой, но без проверки на растворимость. Если же максимальное содержание компонента в первой фазе больше заданной массы, то используется вторая подмодель и определяются только свойства первой фазы по второй подмодели.
Итак, чтобы учитывать эффекты неизотермического процесса, предлагается объединить изотермические модели летучей нефти с жирным газом (PVTO/PVTG) для разных температур в диапазоне, характерном для высокочастотного воздействия, в одну термическую модель. Каждая подобная модель будет в таком случае представлять собой один расчетный шаг неизотермической задачи. Полученное множество термодинамических точек, соответствующих этим моделям, можно расположить в термодинамическом пространстве. PVTO/PVTG подобную модель образуют точки и линейные связи между ними на изотермическом срезе. Если достроить связи между ближайшими точками вне изотермических срезов, образуется единая неизотермическая модель.

Методика расчета свойств фаз, участвующих в смешивающейся фильтрации при ВЧ ЭМ воздействии, состоит в следующем:

• загрузка компонентного состава. Можно использовать практики, являющиеся общепринятыми в нефтяной отрасли;
• настройка на экспериментальные данные для каждой температуры. Очень точная настройка на каждый эксперимент позволяет выбрать необходимую температурную дискретизацию;
• сбор в единую модель с помощью созданного в процессе решения этой задачи python-скрипта для вычисления PVT-свойств в формат симулятора. Расчетный модуль может работать с любой моделью стандартного формата.

Создан python-скрипт для вычисления PVT-свойств на языке программирования Python с использованием библиотек «numpy» и «scipy», где использовался модуль «Interpolation». Получено корректное воспроизведение табличных данных и промежуточных значений для всего диапазона давлений и температуры. В результате получено подтверждение работоспособности данного подхода, который был протестирован на созданной PVT-модели реального месторождения в PVT-симуляторе, для чего были выгружены PVTO/PVTG-модели для температур в диапазоне 50–200 °C, которые были использованы при решении аналитической задачи.