Газовые гидраты представляют собой кристаллические соединения, образующиеся из воды и газа при определенных термобарических условиях. В 40-е годы XX века советскими учеными была выдвинута гипотеза о наличии залежей газовых гидратов в зоне вечной мерзлоты, а в 60-х годах были обнаружены первые месторождения на севере СССР. В настоящее время их продолжают рассматривать в качестве перспективных энергетических ресурсов. В нефтегазодобывающей отрасли нередки случаи образования гидратов в трубопроводах и стволах скважин, в результате чего возникают осложнения и аварийные ситуации. Ликвидация таких последствий гидратообразования требует значительных временных и материальных ресурсов, поэтому основным способом борьбы с нежелательными отложениями гидратов является прогнозирование рисков их возникновения, а также реализация безгидратных режимов работы оборудования на основе точного математического описания условий образования гидратов при моделировании многофазного течения в стволах скважин, трубопроводах, на факельных линиях и т.д. [1–3].

Наличие гидратов значительно увеличивает сложность моделирования фазового равновесия в термодинамической системе (ТДС).
На рисунке 1 показана характерная фазовая диаграмма газоконденсатной системы, содержащей воду. В зависимости от значений давления и температуры в рассматриваемой системе возможно образование одной (V), двух (VL, VA, VH) или трех (VLA, VLH) фаз. На кривой гидратообразования возможно существование и четырех фаз (VLAH V – vapor, L – liquid, A – aqueous, H – hydrate). Стоит отметить, что зона существования гидратов лежит в области положительных температур, а при отрицательных температурах может возникнуть пятая фаза в виде льда.

Оценка состояния равновесия ТДС по кубическим уравнениям состояния (КУС) позволяет описать взаимосвязь между давлением, объемом и температурой многокомпонентных смесей [4]. Развитие КУС привело к разработке более надежных и точных моделей, из которых наиболее часто используемыми оказались уравнения состояния Соаве–Редлиха–Квонга (SRK) и
Пенга–Робинсона (PR). Эти уравнения позволяют с достаточной степенью точности проводить расчеты для оценки фазового состояния многокомпонентных смесей, количества образующихся фаз и их состава. Майкельсен и др. представили систему оценки многофазного равновесия, состоящую из двух этапов: теста на стабильность и расчета фазового равновесия [5]. Данный подход лежит в основе многих коммерческих симуляторов. Несмотря на свою универсальность, расчет многофазного равновесия плохо работает при описании гидратных явлений. Основной причиной является неверное определение летучестей компонент в гидратной фазе, что приводит к ошибочной оценке изменения энергии, необходимой для образования пустых гидратных решеток.

Одной из наиболее широко используемых моделей для прогнозирования равновесия гидратных фаз является модель Ван-дер-Ваальса–Платтеу (VdW–P). Данная модель позволяет определять изменение энергии системы и рассчитывать вероятность заполнения гидратных решеток молекулами газа. Полученные фундаментальные статистические термодинамические уравнения для гидратов определены на основе сходства явления образования гидратов и адсорбции Ленгмюра. В работе [6] Маккой и Синаноглу установили, что модель потенциала Кихары является оптимальным выбором для определения констант Ленгмюра. Однако их расчет довольно таки сложен из-за своей функциональной формы в виде интегральных уравнений. С целью снижения вычислительных затрат Пэриш и Праусниц предложили заменить интегральные уравнения линейными для определения констант Ленгмюра [7].

В данной работе была реализована классическая гидратная модель VdW–P для определения условий образования гидратной фазы с поправками Пэриша–Праусница для определения констант Ленгмюра. Приведен сравнительный анализ полученных результатов с данными стороннего ПО и по корреляции Катца–Карсона [8].

Расчет многофазного равновесия основан на уравнении состояния Пенга–Робинсона с учетом поправки на мольный объем:

где p — давление, R — универсальная газовая постоянная, T — температура, v_e — молярный объем, v — молярный объем с учетом поправки, a, b, b_e, c — параметры уравнения состояния, зависящие от компонентного состава. С помощью уравнения (1), зная давление и температуру, в системе можно определить его мольный объем, а также параметры термодинамической системы: z-фактор, энергию Гиббса, внутреннюю энергию, летучесть компонент и др.

При проведении расчета фазового равновесия выполняется оценка стабильности системы, т.е. расслаивается ли данная смесь на две или более фаз по [5]. Если в системе присутствуют две или более фаз, необходимо, чтобы выполнялось равенство:

где f_i — летучесть i-й компоненты, N_p — число фаз в системе. Это выражение является следствием равновесного состояния термодинамической системы, которая говорит, что летучесть каждой компоненты в различных фазах должна быть одинаковой при равновесии [9].

Кроме выполнения условия (2) должны соблюдаться условия нормировки и материального баланса:

где β₁, β₂, … β_Np — мольные доли фаз, z_i, x_i, x_i,1, x_i,2,… x_i,Np — мольная доля i-й компоненты в исходной и образовавшихся фазах соответственно. С целью уменьшения переменных в уравнениях (2)–(4) вводится коэффициент фазового равновесия, показывающий отношение мольной доли i-й компоненты в одной из фаз к мольной доле той же компоненты контрольной фазы:

где N — число компонент в системе, x_i,F — мольная доля i-й компоненты в контрольной фазе. В качестве контрольной фазы обычно выбирают ту, в которой мольная доля компоненты самая большая (при наличии воды это водная фаза). Решая совместно (3)–(5), можно прийти к следующему соотношению:

где мы получаем N_p – 1 уравнение.

Уравнения (6) являются развитием уравнений Рэчфорда–Райса, решением которых являются количество и доли фаз. Начальные значения для коэффициентов фазового равновесия
определяются с помощью корреляций Вильсона [10], а затем после решения уравнений (6) рассчитываются через коэффициент летучести φ_i,j i-й компоненты в j-й фазе:

Коэффициенты летучести определяются как:

где n_i,j — количество вещества i-го компонента в j-й фазе, V — объем системы, Z — фактор сжимаемости системы. Вычисления проводят до тех пор, пока летучести компоненты во всех фазах не станут равными. Если в системе находится гидратная фаза, то расчет усложняется тем, что выражение (7) не может быть использовано для расчета летучести компонента в гидратной фазе, поскольку уравнение (1) описывает только жидкое и газообразное состояния. Это приводит к возникновению значительных ошибок в расчетах. Кроме того, уравнения (1)–(7) не дают ответа на вопрос об условиях образования гидратов.

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

где μ_w,H — химический потенциал компонента воды в гидратной фазе, μ_w,W — химический потенциал компонента воды в водной фазе. В состоянии равновесия химические потенциалы водной компоненты в фазах гидрата и льда/жидкости равны, т.е.:

Согласно классической модели Ван-дер-Ваальса–Платтеу можно записать:

где индекс 0 относится к тройной точке воды, а индексы H, β и W относятся к гидрату с заполненной решеткой, предполагаемому гидрату с пустой решеткой и воде соответственно, — химический потенциал чистой воды, — разница между химическим потенциалом воды в пустой решетке гидрата и химическим потенциалом жидкой воды в тройной точке, — изменение молярного объема воды между пустой решеткой и жидкой фазой, T̅ — средняя температура между T и T0, — разница в молярной энтальпии воды между пустой решеткой и жидкой фазой, N_hyd — компоненты, участвующие в образовании гидратов, v_k — количество малых (k = 1) и больших(k = 2) полостей, приходящихся на количество молекул воды в одной ячейке, ϴ_ik — парциальная степень заполнения k-й полости i-м компонентом, обозначающая вероятность того, что полость k
будет занята гостевой молекулой i, т.е. молекулы, участвующие в образовании гидратов. Для T₀ и p₀ использованы значения 273,16 К и 611,2 Па соответственно.

Связь химического потенциала чистой воды и воды с растворенным газом можно определить из соотношения [7]:

где f_w,W — летучесть водной компоненты в водной фазе, — летучесть чистой воды. Тогда уравнение (9) с учетом (10) преобразуется к виду:

Значение определяется с помощью следующих уравнений [7]:

где , = = + 0,141(T–T₀) — разница молярной энтальпии и теплоемкости соответственно между чистой жидкой водой и пустой решеткой гидрата в тройной точке.
В таблице 1 приведены основные константы, применяемые в расчете.

Степень заполнения k-й полости i-м компонентом определяется с помощью следующего уравнения:

где C_ik — константы Ленгмюра, полученные с помощью клеточной модели Леннарда–Джонса–Девоншира.

где k_B — постоянная Больцмана, R_k и a — радиусы полости и гидратообразующей молекулы соответственно, w_ik(r) — потенциальная энергия взаимодействия между гостевой молекулой и молекулой воды из полости, r — радиальное расстояние гостевой молекулы от центра полости. Предполагая, что гостевые молекулы имеют сферическую поверхность, энергию взаимодействия можно оценить с помощью потенциала Кихары [12]:

где z — общее количество молекул воды в полости. Значения n, ε, σ, a представляют собой параметры потенциала Кихары, определенные для каждой гостевой молекулы.

Уравнения (14)–(16) представляют собой чрезмерно сложные функциональные зависимости. Для снижения вычислительных ресурсов Пэрриш и Праусниц предложили заменить интегральное уравнения (14) следующим выражением [7]:

где A_ik, B_ik — постоянные, определяемые из табличных данных.

Оценку летучестей гостевых компонент в гидратной фазе Майкельсон предложил определять из внутренней итерационной схемы, используя соотношения между степенью заполнения полостей гостями-молекулами и общим количеством пустых полостей. Для этого необходимо решить уравнение [10]:

где Y — отношение количества пустых больших полостей к общему количеству пустых полостей в ячейке гидрата, — соотношение между мольной долей молекул типа i и мольной долей воды в гидратной структуре, — общее количество пустых полостей на моль воды. Тогда летучесть гостевых молекул равна:

Соотношения (1), (2), (5)–(8), (11)–(13), (17)–(19) позволяют получить условия образования гидратов по заданному компонентному составу и давлению в равновесной ТДС. Стоит учитывать, что данный процесс является итерационным, т.к. изначально при переходе воды в гидрат летучести компонент в гидратной фазе не известны. Постепенно в процессе последовательных вычислений, изменяя температуру, получаем равенство летучестей во всех фазах для заданного давления.

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

с точностью до ε = 10^-13.

Разработанный алгоритм был внедрен в корпоративные ПО «РН-СИМТЕП» [13]
и ПК «РН-ВЕГА» [14–16]. Результаты расчетов сравнивались с данными стороннего ПО [17], экспериментальными данными и значениями, полученными по корреляции Катца–Карсона.
На рисунке 2 показаны результаты сравнения данных по чистым компонентам и смеси газа. Как видим, представленная модель достаточно хорошо описывает условия гидратообразования, а результаты согласуются с данными стороннего ПО. В таблице 2 указаны средние отклонения значений давлений, полученных с помощью различных ПО, от давлений, полученных в ходе эксперимента.

Для чистых газов результаты эксперимента были взяты из работы [18]. Экспериментальные данные по гидратообразованию для промыслового газа, состав которого представлен в таблице 3, получены в результате лабораторных исследований.

Как и предполагалось, ячеистая модель дает результаты с погрешностью в пределах 10 %. Авторы статьи [19] утверждают, что при обработке более 1500 наборов данных они также заметили, что рассматриваемая модель для чистых газов работает хорошо, однако для смесей погрешность увеличивается в 2–3 раза. Кроме того, не стоит исключать и погрешность измерения экспериментальных данных по гидратообразованию.