Расчет фазовых равновесий чистых легких углеводородов с помощью уравнения состояния PC-SAFT и алгоритма прямой минимизации энергии Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

ФИЗИКА ЗЕМЛИ, АТМОСФЕРЫ И ГИДРОСФЕРЫ

Расчет фазовых равновесий чистых легких углеводородов с помощью уравнения состояния PC-SAFT и алгоритма прямой минимизации энергии

А. В. Исаева,1' * М.О. Бевзо1'f

1 Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра физики Земли Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2 (Поступила в редакцию 04.10.2023; после доработки 20.10.2023; подписана в печать 23.10.2023)

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

PACS: 64.30.+t, 64.70.Fx. УДК: 536.71, 536.77, 544.344.

Ключевые слова: алгоритм прямой минимизации энергии, уравнение состояния PC-SAFT, уравнение состояния Пенга—Робинсона, задача линейного программирования.

DOI: 10.55959/MSU0579-9392.78.2360903

введение

Прогноз фазовых превращений природных неф-тей и газоконденсатных смесей необходим для решения широкого круга инженерных задач нефтегазовой отрасли. Для расчета параметров фаз используют различные уравнения состояния [1]. Наиболее популярны кубические уравнения состояния: уравнение Пенга-Робинсона, уравнение Со-аве-Редлиха-Квонга, уравнение Зудкевича-Иоффе и др. Перечисленные уравнения состояния доступны в промышленных программных продуктах — композиционных симуляторах [2, 3]. Композиционные симуляторы предназначены для расчета течений углеводородных смесей в залежах с учетом возможных фазовых превращений и межфазного массообмена. К примеру, уравнение состояния Пенга-Робинсона имеет вид

Р

RT

v — b v(v + b) + b(v — b)'

(1)

где p — давление, v — мольный объем фазы, T — температура, R — универсальная газовая постоянная, параметры a и b зависят от критических температур, критических давлений и других параметров рассматриваемых веществ или их смесей [1].

В то же время растет интерес к использованию уравнения состояния PC-SAFT для описания фазового поведения углеводородов. В работах [4] и [5] показано, что это уравнение состояния дает преимущества в точности прогноза плотности и изобарной теплоемкости природных углеводородных смесей

* E-mail: avisaeva@mail.ru t E-mail: bevzomo@gmail.com

по сравнению с кубическими уравнениями состояния. Уравнение состояния PC-SAFT задает выражение для свободной энергии Гельмгольца рассматриваемой системы (вещества) [6]. При этом энергия Гельмгольца F определена суммой вкладов энергий различной природы:

F Fhc + Fpert?

где Fhc — энергия в системе твердых цепей, Fpert — поправка к этой энергии, рассчитываемая по теории возмущений [6]. Для чистого вещества энергия Fhc определяется выражением:

Fhc = Fid + mFhs — (m — 1) ln ghs(a),

где Fid — энергия Гельмгольца одного моля одно-компонентного идеального газа, Fhs — энергия в системе твердых сфер, #hs(a) — радиальная функция распределения для системы твердых сфер, зависящая от диаметра сегмента молекулы а рассматриваемого вещества, m — число сегментов в молекуле.

Энергию в системе твердых сфер можно представить в виде

Fhs =

+

Z3

ЗС1С2

1 — Сэ Сз (1 - Сз):

Z3

г +

+ [ Ц -Со J 1X1(1 -Сз)

здесь Сп = прт!п/ 6, п € {0,1, 2, 3}, где р = Жд/V, Жд — число Авогадро, ! — диаметр сегмента молекулы с учетом его температурной зависимости:

! = а [1 - 0.12 ехр (—3 е/ кТ)],

е — энергетический параметр сегмента, к — постоянная Больцмана.

а

Радиальная функция распределения выражается как

ghs(^)

1

+

3dZ2

d2 Z2

+

1 - Zs 2(1 - Сз)2 2(1 - Сз)

з •

В рассчитываемую по теории возмущений поправку -Рре!^ к энергии в системе твердых цепей входит вклад обусловленный действием дис-

персионных сил. Этот вклад может быть записан в виде

^¡вр = -2тгр11(1],т.)т2 ) сг3-

- -кртС-хЬ^Ъ т)т2

где п = Сз, 11(п,т) и 12(п,т) — интегралы возмущений, С1 — сжимаемость.

Интегралы возмущений можно приближенно записать в виде степенных рядов по п:

6

Ii(n, m) = ^ ai(m)rf

i=0 6

I2(n,m) = bi(m)rf

где а^ и Ъ,, — коэффициенты, которые можно выразить в виде, предложенном в статье [7]:

ai (m) = aoi + bi(m) = boi +

m - 1

m m-1

m

aii + bii +

(m - 1)(m - 2)

m2

(m — 1)(m — 2)

a2i,

b2i-

Коэффициенты а^, a1i, a2i, Ъ0.1, Ъ1i, Ъ^ определяют экспериментально. Конкретные значения параметров уравнений состояния, используемые далее в расчетах, представлены в Приложении к настоящей статье.

Сжимаемость С1 задается выражением

Ci

1 + m

8n - 2n2

+

+(1 -

(1 - n)4

20n - 27n2 + 12ns - 2n4"

[(1 - n) (2 - n)]2

Подробнее об уравнении состояния PC-SAFT можно узнать из литературы [4]—[7].

Заметим, что энергия Гельмгольца изохорно-изо-термической системы (заключенной в фиксированном объеме при постоянной температуре) достигает минимума в состоянии термодинамического равновесия [1]. Этот факт используется в методе прямой минимизации энергии для расчета параметров равновесного состояния. Как будет показано далее, в дискретной постановке нелинейная задача минимизации энергии системы сводится к задаче линейного программирования. Такой подход хорошо известен в геохимии и петрологии, где он активно

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

В настоящее время метод прямой минимизации энергии находит применение в нефтегазовой отрасли. С его помощью рассчитывают параметры парожидкостного равновесия углеводородных смесей [11], а также трехфазных равновесий смесей углеводородов и диоксида углерода [12]. Поскольку уравнение состояния PC-SAFT в явном виде определяет энергию Гельмгольца вещества или смеси веществ, представляется перспективным использовать это уравнение состояния в сочетании с прямой минимизацией энергии Гельмгольца для расчета параметров фазовых равновесий. В настоящей работе мы демонстрируем продуктивность такого подхода для расчета параметров парожидкостно-го равновесия легких углеводородов (метана, этана, пропана и н-бутана). Для легких углеводородов параметры парожидкостного равновесия (давления насыщения, мольные объемы равновесных фаз) определены экспериментально для ряда температур [13]-[16]. Это позволяет сопоставить расчетные значения параметров парожидкостного равновесия со значениями, известными из лабораторных экспериментов.

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

1. алгоритм прямой минимизации энергии

Рассмотрим один моль чистого вещества. Согласно правилу фаз Гиббса [1], в состоянии термодинамического равновесия чистое вещество может образовывать не более трех фаз. Обозначим число фаз Ыр, следовательно, Ыр может принимать значения 1, 2 или 3. Если рассматривать только парожид-костные равновесия, то максимальное число фаз сокращается до двух.

2

m

i

Если доступный системе объем гнун и температура Т остаются постоянными, то в состоянии термодинамического равновесия свободная энергия Гельмгольца системы достигает минимума [1]. Запишем это условие математически для многофазной системы:

N

^ =

N

шт

Е^а (ТО

(2)

где Еа (Т, га) — мольная энергия Гельмгольца фазы а, которая зависит от температуры Т и мольного объема фазы га, па — молярная доля этой фазы. Поскольку мы рассматриваем один моль вещества, будут справедливы условия

N

У^Пд = 1, 0 < Па < 1.

(3)

Условие постоянства общего объема гнун многофазной системы примет вид

N

Е

а=1

^аПа — ^вув.

(4)

Тогда параметры равновесного состояния системы (количество фаз, их доли и мольные объемы) можно рассматривать как решение задачи минимизации (2) с учетом ограничений (3) и (4):

объема фаз из всех г^, г — 1,..., Ы^, а также соответствующие им ненулевые значения долей фаз па. Тогда задачу (6) можно переписать эквивалентно как

/ N \а=1

(7)

где п — набор значений, выражающих доли допустимых фаз, т.е. мольные объемы этих фаз содержатся в наборе г^, г — 1,...,Ы^. Следовательно, в результирующем наборе значений п*, который будет соответствовать минимуму энергии Гельм-гольца системы (7), ненулевыми будут одно или два значения, определяющие молярные доли фаз в равновесии.

Заметим, что сумму произведений мольных энергий Гельмгольца фаз с их долями в (7) можно воспринимать как скалярное произведение пары векторов в Ы^-мерном пространстве. Поэтому задачу о минимизации скалярного произведения (7) с ограничениями (3) и (4) можно рассматривать как задачу линейного программирования [17].

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

N

Ыр,г*а,пд — а^шт Е ^а(Т,

, (5)

где Ы* — количество фаз в равновесии, гд — мольный объем равновесной фазы а, пд — ее молярная доля. Таким образом, получаем, вообще говоря, нелинейную задачу минимизации с ограничениями.

Для поиска численных решений задачи (2)—(4) можно воспользоваться следующим приемом. Обозначим гт1п и гтах наибольшее и наименьшее допустимое значение мольного объема фазы соответственно. К примеру, при использовании уравнения состояния Пенга-Робинсона (1) естественным ограничением снизу для величины мольного объема фазы является параметр 6, т.е. можем считать гтш — 6.

Ограничимся далее случаем парожидкостных систем (Ыр может принимать значения 1 или 2). Зададим сетку на множестве допустимых значений мольного объема фазы: г^ € (гт;п,гтах), г — 1,...,Ы^, где Ы — количество узлов сетки. Тогда условие минимальности свободной энергии Гельмгольца системы (2) можно записать как

N

шт

N

Е^а(ТО

(6)

1

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

Рис. 1. Блок-схема алгоритма расчета параметров фазового равновесия чистого вещества (УС — уравнение состояния)

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

а

а

и™ )П

а

0.8 г

0 5 10 15

Приведенный мольный объем

Рис. 2. Изотерма уравнения состояния Пенга-Робинсо-на: 1 — зависимость приведенного давления от приведенного объема; 2 — физически нереализуемый участок изотермы; 3 — давление парожидкостного равновесия)

6

0 5 10 15

Приведенньт мольный объем

Рис. 3. Зависимость приведенной энергии Гельмголь-ца от приведенного мольного объема: 1 — энергии Гельмгольца в соответствии с уравнением состояния Пенга-Робинсона; 2 — физически нереализуемый участок; 3 — энергия Гельмгольца двухфазной системы

Приведенный мольный объем системы

■е

г

и (в ю о

15

10

5

0 5 10 15

Приведенный мольный объем системы

Рис. 4. Результат расчета параметров равновесных фаз: а — зависимость молярных долей фаз от объема системы; б — молярные объемы фаз; 1 — жидкая фаза; 2 — газовая (паровая) фаза

ют в равновесии при рассматриваемых температуре и объеме системы.

Рис. 2-5 иллюстрируют проводимые расчеты на примере метана и уравнения состояния Пенга-Робинсона. На рис. 2 показана изотерма этого уравнения состояния: зависимость приведенного (безразмерного) давления от приведенного мольного объема при фиксированной температуре. Давление отнесено к величине критического давления для метана, а мольный объем — к параметру Ь уравнения состояния Пенга-Робинсона. Рассматривалась температура ниже критической температуры для метана. Штриховой линией показан участок изотермы, физически нереализуемый для чистого вещества [1].

На рис. 3 для той же температуры показана зависимость от приведенного мольного объема сво-

бодной энергии Гельмгольца, отнесенной к произведению температуры на универсальную газовую постоянную. Заметим, что между точками, отмеченными ромбами на рис. 3, кривая энергии Гельм-гольца идет выше прямой, соединяющей эти две точки. Физическая интерпретация описанного факта следующая: в этой области мы наблюдаем двухфазное равновесное состояние системы, поскольку энергия Гельмгольца двухфазной системы оказывается меньше энергии однофазной. Действительно, энергии двухфазной системы будут соответствовать точки на прямой, т.к. ее энергия складывается из произведений мольных энергий Гельмгольца фаз с их долями (2). Эта энергия оказывается меньше энергии Гельмгольца однофазной системы, задаваемой уравнением состояния. Следовательно, точки, отмеченные ромбами на рис. 3, соответствуют реше-

Я 0.7

0.6

0.5

Я

0.4

0 5 10 15

Приведенный мольный объем системы

£7.6

0 5 10 15

Приведенный мольный объем системы

Рис. 5. Результат расчета параметров равновесной системы: а зависимость энергии Гельмгольца системы от ее объема

зависимость давления от объема системы; б

нию задачи минимизации (7) с ограничениями (3) и (4): проекции этих точек на ось абсцисс определяют мольные объемы равновесных фаз. Мольный объем жидкой фазы меньше мольного объема газовой (паровой) фазы. Заметим, что для таких фаз давление одинаково и соответствует давлению паро-жидкостного равновесия, что иллюстрирует рис. 2.

Таким образом, для определения параметров па-рожидкостного равновесия мы поступали следующим образом. При фиксированной температуре увеличивали объем системы vsys начиная с vmin = b. Когда vsys становится больше мольного объема равновесной жидкой фазы, в решении задачи линейного программирования (3), (4), (7) появляются два ненулевых значения для молярных долей равновесных фаз. С дальнейшим ростом vsys молярная доля газовой фазы растет, а жидкой — уменьшается. При превышении vsys мольного объема равновесной паровой фазы жидкая фаза исчезает. На рис. 4, а представлены такие молярные доли фаз, получаемые в результате решение задачи линейного программирования, в зависимости от vsys. Рис. 4, б показывает зависимость молярного объема фаз от vsys. Ясно, что в однофазных областях найденный в качестве решения задачи линейного программирования мольный объем фазы совпадает с vsys. При этом в двухфазной области получим два решения — мольные объемы жидкой и газовой фаз.

Рис. 5, а показывает найденную зависимость давления от vsys, а рис. 5, б — зависимость свободной энергии Гельмгольца от vsys. Эти зависимости соответствуют представленным на рис. 2 и 3. Как обсуждалось выше, всюду в двухфазной области мы наблюдаем одно и то же давление парожидкостного равновесия.

Сделаем два важных замечания. Во-первых, сведение исходной задачи минимизации (2)—(4) к задаче линейного программирования (3), (4), (7) существенно упрощает программную реализацию обсуждаемого алгоритма. В популярных языках программирования доступны готовые функции, поз-

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

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

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

Параметры парожидкостного равновесия рассчитывались для четырех легких алканов: метана, этана, пропана и н-бутана. Каждое вещество рассматривалось при шести различных температурах. Для рассмотренных температур опубликованы давления насыщения (парожидкостного равновесия) и мольные объемы фаз, известные для этих веществ из лабораторных экспериментов [ 13]—[16]. Таким образом, удалось сопоставить опубликованные параметры с рассчитанными в настоящей работе.

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

Таблица 1. Параметры парожидкостного равновесия метана при разных температурах

Температура, К Мольный объем жидкой фазы, Мольный объем газовой фазы, Давление парожидкостного

мл/моль мл/моль равновесия, МПа

42.565 ± 0.009 (эксп.) 1578.9 ± 0.4 (эксп.) 0.64283 ± 0.00011 (эксп.)

140 39 (ПР) 1582 (ПР) 0.641 (ПР)

43 (PC-SAFT) 1605 (PC-SAFT) 0.640 (PC-SAFT)

47.697 ± 0.010 (эксп.) 631.62 ± 0.15 (эксп.) 1.5929 ± 0.0002 (эксп.)

160 45 (ПР) 622 (ПР) 1.600 (ПР)

47 (PC-SAFT) 646 (PC-SAFT) 1.586 (PC-SAFT)

58.078 ± 0.012 (эксп.) 261.14 ± 0.07 (эксп.) 3.2868 ± 0.0004 (эксп.)

180 59 (ПР) 252 (ПР) 3.306 (ПР)

59 (PC-SAFT) 267 (PC-SAFT) 3.272 (PC-SAFT)

В табл. 1 такое сопоставление представлено на примере метана для трех температур. В этой таблице отметка «(эксп.)» около значения указывает, что оно получено в лабораторном эксперименте [13]. Отметки «(ПР)» и «(PC-SAFT)» обозначают, что при расчете параметра использовалось уравнение состояния Пенга-Робинсона или PC-SAFT соответственно. Все расчеты проводились с использованием разработанного авторами статьи комплекса программ в GNU Octave. При расчетах параметры уравнения состояния Пенга-Робинсона, а также критические температуры, давления, ацентрические факторы чистых веществ были взяты из работы [18]. Для расчетов с использованием уравнения состояния PC-SAFT параметры легких алканов (число сегментов в молекуле, диаметр сегмента молекулы, энергетический параметр сегмента молекулы) были взяты из статьи [19]. При этом нами использовалась матрица универсальных коэффициентов, необходимых в уравнении состояния PC-SAFT для расчета интегралов возмущений, из работы [4]. Конкретные значения параметров уравнений состояния и углеводородов, использованные в расчетах, представлены в Приложении к настоящей статье.

Из табл. 1 видно, что погрешности экспериментальных данных пренебрежимо малы в сравнении с той разницей, которая получается при сопоставлении расчетных значений и экспериментальных данных. Действительно, относительное отклонение расчетных и экспериментальных значений достигает 8%, в то время как относительная погрешность экспериментальных данных оказывается менее десятых долей процента. Поэтому мы можем беспрепятственно принять экспериментальные данные за «эталон» при сравнении результатов расчетов параметров парожидкостного равновесия, полученных с использованием разных уравнений состояния. Это помогает оценить, с какой точностью уравнения состояния предсказывают фазовое поведение углеводородов.

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

в прикладных задачах (200-500 К). Это обстоятельство продиктовано целью сопоставить рассчитанные нами параметры парожидкостного равновесия легких алканов со значениями, известными из лабораторных экспериментов [13]—[16]. Соответственно, мы рассматриваем вещества при тех же температурах, при которых они изучались в состоянии парожидкостного равновесия в работах [13]-[16]. Разумеется, представляемые нами способы расчета фазовых равновесий можно использовать и в других температурных диапазонах в зависимости от решаемой задачи.

Для удобства сравнения результатов, полученных для различных веществ, рассматривались относительные отклонения (в процентах) расчетных параметров от их экспериментальных значений из литературы [13]-[16]. Изучались относительные отклонения расчетных и экспериментальных значений мольного объема жидкой фазы, мольного объема газовой (паровой) фазы и давления парожид-костного равновесия.

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

Рис. 6-8 демонстрируют результаты расчетов, полученные для метана. На рис. 6 показана зависимость от приведенной температуры относительных отклонений расчетных и экспериментальных значений для мольного объема жидкой фазы, на рис. 7 — для мольного объема газовой фазы, на рис. 8 — для давления парожидкостного равновесия. Закрашенные черным цветом маркеры отражают результаты расчетов, полученные с использованием уравнения состояния Пенга-Робинсона, светлые маркеры — результаты, полученные с использованием уравнения состояния PC-SAFT. Видно, что для метана оба уравнения состояния демонстрируют хорошую точность прогноза параметров парожидкост-ного равновесия: относительные отклонения рас-

О Уравнение PC-SAFT • Уравнение Пенга-Робинсона

О

J_L_nl_

0,65 0,70 0,75 0,80 0,85 0,90 0,95 Приведенная температура

о ">?

2 .

а з

и се

0Í 3

§ «

§ я

о

О Уравнение PC-SAFT • Уравнение Пенга-Робинсона

О О

0,65 0,70 0,75 0,80 0,85 0,90 0,95 Приведенная температура

Рис. 6. Относительные отклонения расчетных значений мольного объема жидкой фазы метана от соответствующих экспериментальных значений в зависимости от приведенной температуры (УС ПР — уравнение состояния Пенга—Робинсона, УС PC-SAFT — уравнение состояния PC-SAFT)

Рис. 7. Относительные отклонения расчетных значений мольного объема газовой фазы метана от соответствующих экспериментальных значений в зависимости от приведенной температуры (УС ПР — уравнение состояния Пенга—Робинсона, УС PC-SAFT — уравнение состояния PC-SAFT)

О Уравнение PC-SAFT • Уравнение Пенга-Робинсона

О

_ • О

° 0 о

J_I_I_I_I_I_I_I_I_I_I

' 0,65 0,70 0,75 0,80 0,85 0,90 0,95 Приведенная температура

Рис. 8. Относительные отклонения расчетных значений давления парожидкостного равновесия метана от соответствующих экспериментальных значений в зависимости от приведенной температуры (УС ПР — уравнение состояния Пенга—Робинсона, УС PC-SAFT — уравнение состояния PC-SAFT)

четных и экспериментальных значений не превышают 11%. Заметим, что в расчетах с использованием уравнения состояния Пенга-Робинсона относительные отклонения для мольных объемов жидкой фазы получились в среднем больше, чем относительные отклонения для мольных объемов газовой фазы. Такой результат для этого уравнения состояния

ожидаем и описан в [1].

На рис. 9-11 показаны аналогичные результаты расчетов для четырех легких углеводородов парафиновой группы: метана, этана, пропана и бутана. Видно, что для всех этих веществ рассмотренные уравнения состояния демонстрируют сопоставимую точность прогноза параметров парожид-костного равновесия. Относительные отклонения мольного объема жидкой фазы, полученные при использовании уравнения состояния PC-SAFT, не превышают 6%, а при использовании уравнения состояния Пенга-Робинсона — 11%. Таким образом, уравнение состояния PC-SAFT точнее предсказывает величины мольных объемов жидкой фазы. В то же время уравнение состояния PC-SAFT демонстрирует меньшую точность прогноза мольных объемов газовой фазы (наблюдаются относительные отклонения до 16%) и давления парожидкостного равновесия (относительные отклонения — до 5%), чем уравнение состояния Пенга-Робинсона.

Заметим, что на рис. 9-11 для уравнения состояния PC-SAFT прослеживаются некоторые общие тенденции. Во-первых, для веществ с большей молекулярной массой относительные отклонения расчетных значений параметров парожидкост-ного равновесия от их экспериментальных значений выше (наблюдаемые отклонения для пропана и бутана в среднем выше, чем для метана и этана). Во-вторых, для всех рассмотренных веществ прогноз значений мольных объемов газовой фазы и давления парожидкостного равновесия ухудшается с ростом температуры (приближение к критической температуре вещества). При этом подобная тенденция не прослеживается для уравнения

1.2

§ 1,0

& и

¡ц

а §

I и

0,8 о. о

8 0,6

0,4

о &

в 0,2

16 14 12 ,10

cd 6

S °

<L>

i 4

Метан (ПР) Метан (PC-SAFT) Этан (ПР) Этан (PC-SAFT) Пропан (ПР) Пропан (PC-SAFT) н-Бутан (ПР) н-Буган (PC-SAFT)

is ▲

О

_L

_L

z

i

JñL

0,65 0,70 0,75 0,80 0,85 0,90 0,95 1,00 Приведенная температура

s*

Si я

§ 1 s ^ g >1 H 1

о 3

■ Метан (ПР)

□ Метан (PC-SAFT)

• Этан (ПР)

с- Этан (PC-SAFT)

А Пропан (ПР) J

. л Пропан (PC-SAFT)

★ н-Буган (ПР)

А н-Буган (PC-SAFT) * д

д °

О О

а Л

_L

_L

l'

_L

_L

J

0,65 0,70 0,75 0,80 0,85 0,90 0,95 1,00 Приведенная температура

Рис. 9. Относительные отклонения расчетных значений мольного объема жидкой фазы углеводородов от соответствующих экспериментальных значений в зависимости от приведенной температуры (ПР — уравнение состояния Пенга—Робинсона, PC-SAFT — уравнение состояния PC-SAFT)

Рис. 10. Относительные отклонения расчетных значений мольного объема газовой фазы углеводородов от соответствующих экспериментальных значений в зависимости от приведенной температуры (ПР — уравнение состояния Пенга—Робинсона, PC-SAFT — уравнение состояния PC-SAFT)

« 3

о. о

Метан (ПР) Метан (PC-SAFT) а Этан (ПР) Этан (PC-SAFT)* Пропан (ПР) Пропан (PC-SAFT) н-Буган (ПР) н-Буган (PC-SAFT)

¿t *

A Í3

-О О

'et □

.*!»» I«

0,60 0,65 0,70 0,75 0,80 0,85 0,90 0,95 1,00 Приведенная температура

Рис. 11. Относительные отклонения расчетных значений давления парожидкостного равновесия углеводородов от соответствующих экспериментальных значений в зависимости от приведенной температуры (ПР — уравнение состояния Пенга—Робинсона, PC-SAFT — уравнение состояния PC-SAFT)

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

SAFT имеет другой способ параметризации — чистое вещество определяется заданием числа сегментов в молекуле, диаметром сегмента молекулы и его энергетическим параметром.

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

ЗАКЛЮЧЕНИЕ

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

тических вычислений GNU Octave, в которой программно реализовали описанные в настоящей статье алгоритмы. Созданные программы позволяют моделировать фазовое поведение чистых веществ с использованием уравнения состояния PC-SAFT, а также с использованием одного из наиболее востребованных в нефтегазовой отрасли уравнений состояния — уравнения Пенга-Робинсона. Поэтому с помощью разработанных программ удалось сопоставить точность прогнозирования параметров па-рожидкостного равновесия легких углеводородов с использованием этих двух уравнений состояния.

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

личных температурах. Вычислялись относительные отклонения расчетных и экспериментальных значений параметров парожидкостного равновесия (мольных объемов жидкой и газовой фазы, а также давления парожидкостного равновесия). Для рассмотренных веществ и температур получена хорошая точность прогноза параметров парожидкостного равновесия — относительные отклонения не превышают 16%. При этом уравнение состояния PC-SAFT оказалось точнее в прогнозировании значений мольных объемов жидкой фазы, а уравнение состояния Пенга-Робинсона — в прогнозировании значений мольных объемов газовой фазы и давления парожидкостного равновесия. Таким образом, в настоящей статье мы продемонстрировали возможность использования метода прямой минимизации энергии в сочетании с уравнением состояния PC-SAFT для расчета параметров парожидкостного равновесия легких углеводородов.

приложение

Таблица 2. Параметры легких углеводородов [18, 19]

Вещество Тс, К Рс, МП а со m а, А e/fe, К

Метан 190.56 4.599 0.0115 1.0000 3.7039 150.03

Этан 305.32 4.872 0.0995 1.6069 3.5206 191.42

Пропан 369.83 4.248 0.1523 2.0020 3.6184 208.11

н-Бутан 425.12 3.796 0.2002 2.3316 3.7086 222.88

Таблица 3. Матрица универсальных параметров уравнения состояния PC-SAFT [4]

i aoi ai i «2 ¿

0 0.876214348719 -0.356138072415 -0.033972894794

1 1.359414877196 1.197006547003 -1.027538490455

2 -3.753331119981 -11.077288420480 12.221043599616

3 3.200447627334 60.337470085961 -48.551834976584

4 21.162143755838 -173.091309439975 108.459219986741

5 -53.542254106988 261.362994582518 -140.715575522984

6 29.331248242432 -164.459042961042 81.499291851005

i boi ^lz &2¿

0 0.679333892140 0,597244160061 0.038554880712

1 3.063826512650 1.237865026919 1.061542913488

2 -12.625108297584 -8.331327739620 0.424275193059

3 21.047164840594 61.456368121717 -56.307434168993

4 -45.402862229694 -44.972742056517 143.087511514563

5 180.515819639032 122.262049217030 -19.703093779073

6 -231.492587269449 -196.635009485523 -146.631916746802

В уравнении состояния Пенга-Робинсона (1) па- где Д — универсальная газовая постоянная, Тс — раметры а и Ь рассчитывались как критическая температура вещества, Рс — его кри-

тическое давление, параметр О а рассчитывался по

а = О а ■ (Д ■ Тс)2/Рс, Ь = Ов0 (Д ■ Тс)/Рс,

формуле

ПА = ПАо (1 + (0.37464+ 1.54226 ш - 0.26992 ш2) х

х(1 - (Т/ТС)0'5))2,

в которой ш — ацентрический фактор вещества, ПАо = 0.457235, ПБо = 0.077796.

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

В уравнении состояния PC-SAFT вещество параметризуется числом сегментов в молекуле, диаметром сегментов молекулы и энергетическим параметром сегментов молекулы. Конкретные значения этих параметров, использованные в расчетах, также представлены в табл. 2.

Как обсуждалось во Введении, в уравнении состояния PC-SAFT расчет интегралов возмущений зависит от выбора значений коэффициентов (универсальных параметров) aoi, ац, a2i, boi, Ъц, &2г. В настоящей работе для расчетов использовались значения этих коэффициентов, предложенные в работе [4]. Конкретные значения коэффициентов представлены в табл. 3.

[2

[3

[4 [5 [6 [7 [8 [9 [10 [11

Брусиловский А.И. // Фазовые превращения при разработке месторождений нефти и газа. М.: Грааль, 2002.

Schlumberger Limited. ECLIPSE Technical Description. Houston, TX, USA: Schlumberger Limited, 2014.

Schlumberger Limited. INTERSECT Reservoir Simulator Technical Description. Houston, TX, USA: Schlumberger Limited, 2012.

Александров И.С., Григорьев Б.А. // Вести газовой науки. № 4 (36). 237 (2018).

Григорьев Б.А., Александров И.С., Герасимов А.А. // Газовая промышленность. № 6 (769). 52. (2018). Gross, J., Sadowski, G. // Fluid Phase Equilibria. 2(168). 183. (2000).

Liu H, Hu. Y. // Fluid Phase Equilibria. 122(1-2). 75. (1996).

Connolly J.A.D. // Earth and Planetary Science Letters. 236(1-2), 524 (2005).

Коннолли Дж.А.Д. // Петрология. 25(5). 533. (2017).

Vrijmoed J.C., Podladchikov Y.Y. // Geochemistry,

Geophysics, Geosystems. 23, 4. (2022).

Исаева А.В., Доброжанский В.А., Хакимова Л.А.,

[12]

13

14

15

16

17

18

19

Подладчиков Ю.Ю. // Газовая промышленность. № 2 (812). 22. (2021).

Khakimova L., Isaeva A., Dobrozhanskiy V., Podladchikov Y. Direct Energy Minimization Algorithm for Numerical Simulation of Carbon Dioxide Injection / SPE Russian Petroleum Technology Conference, Virtual, October 2021. SPE-206611-MS.

Kleinrahm R., Wagner W. // The Journal of Chemical Thermodynamics. 8(18). 739. (1986). Douslin D.R., Harrison R.H. // The Journal of Chemical Thermodynamics. 5(4). 491. (1973). Thomas R.H.P., Harrison R.H. // Journal of Chemical & Engineering Data. 27(1). 1. (1982). Kay W.B. // Industrial & Engineering Chemistry. 32(3). 358. (1940).

Васильев Ф.П. Методы оптимизации. М.: Изд. «Факториал Пресс», 2022.

Danesh A. PVT and Phase Behaviour Of Petroleum Reservoir Fluids // Elsevier Science B.V. 1998. 47. 1-388. eBook ISBN: 9780080540054 Gross J., Sadowski G. // Ind. Eng. Chem. Res. 40. 1244. (2001).

1

Calculation of Phase Equilibria of Pure Light Hydrocarbons Using the PC-SAFT Equation of State and the Direct Energy Minimization Algorithm

A. V. Isaeva1", M. O. Bevzo16

1 Faculty of Physics, Lomonosov Moscow State University. Moscow 119991, Russia E-mail: aavisaeva@mail.ru, bbevzomo@gmail.com

In the oil and gas industry, there is growing interest in using the PC-SAFT equation of state for predicting the phase behaviour and physical properties of hydrocarbons. This article demonstrates how this equation of state can be used in combination with the direct energy minimization algorithm of an isochoric-isothermal system for calculating vapour-liquid equilibrium parameters of hydrocarbons. The proposed approach is tested on four substances: methane, ethane, propane, and n-butane.

PACS: 64.30.+t, 64.70.Fx.

Keywords: direct energy minimization algorithm, PC-SAFT equation of state, Peng-Robinson equation of state, linear programming problem. Received 04 October 2023.

English version: Moscow University Physics Bulletin. 2023. 78, No. 6. Pp. 876-886. Сведения об авторах

1. Исаева Анна Вячеславновна — канд. физ.-мат. наук, доцент; e-mail: avisaeva@mail.ru.

2. Бевзо Мария Олеговна — студент; e-mail: bevzomo@gmail.com.