Теплофизика высоких температур, 2021, T. 59, № 6, стр. 852-859

ВВЕДЕНИЕ

Интерес к исследованию поведения веществ в около- или закритических состояниях (сверхкритические флюиды – СКФ), представляющих нечто среднее между жидкостью и газом, вызван их специфическими свойствами. Среди них важно отметить лишь то, что они могут сжиматься как газы (обычные жидкости практически несжимаемы) и в то же время способны растворять твердые вещества, что газам не свойственно [1]. Сама критическая точка (КТ), определяющая границу СКФ, является фундаментальным физическим параметром вещества, не менее важным, чем точки плавления или кипения [2]. Поэтому измерения и расчеты как термодинамических свойств вещества в СКФ, так и критических параметров являются одной из основных задач теплофизики с конца 19-го века, начиная со знаменитого уравнения Ван-дер-Ваальса [3]. И к настоящему времени для многих газов и жидкостей на основе этих исследований были созданы базы данных и построены так называемые референтные уравнения состояния (РУРС), например, [4]. Очевидно, что закритическое состояние и КТ находились в сфере интересов В.Е. Фортова на протяжении практически всей его научной деятельности. Свидетельством этому является обзор по координатам КТ, написанный им с его коллегами более 40 лет назад [5], а также многочисленные монографии и обзоры по уравнениям состояния самых различных систем в широком диапазоне параметров, которые создавались ими до самого недавнего времени (см., например, [6–9]). Несмотря на прогресс в этой области, стоит отметить, что даже РУРС носят сугубо эмпирический характер. Поэтому область их применения по плотностям, температурам, давлениям и т.д. ограничена только тем диапазоном, где имеются надежные измерения, данные которых и используются при расчете эмпирических констант. Для расчетов вне этого диапазона необходимо экстраполировать данные, полученные по этим эмпирическим уравнениям, а такие экстраполяции могут приводить к серьезным ошибкам, если уравнение состояния не имеет теоретической основы. Кроме этого, можно использовать различные соотношения подобия [10, 11]. Такой подход имел место в обзоре [5]. Соотношения подобия очень удобны, но при построении уравнения состояния (УРС) их тоже желательно объединять с более строгими теоретическими методами. Именно в этом и состоит цель настоящей работы.

Одним из строго обоснованных и общих теоретических методов для построения УРС является вириальное разложение. Оно абсолютно точно в области сходимости этого разложения [2]. Более того, возможно построить и аналитическое продолжение вириального ряда за пределы этой области [12]. Но на пути этого подхода встают серьезные трудности, связанные с вычислением старших вириальных коэффициентов, их нерегулярным поведением, а также плохой сходимостью вириального ряда. Тем не менее за последнее время наметился определенный прогресс при вычислении вириальных коэффициентов высокого порядка. Так, например, для системы частиц, описываемых потенциалом Леннарда−Джонса (ЛД), в настоящее время вычислено 16 коэффициентов [13]. Но погрешность старших (начиная с 10-го) все еще высока. Обычно стандартный вириальный ряд для давления строится в виде разложения по степеням стремящейся к нулю плотности и первый его член дает давление идеального газа. Однако при повышении плотности для описания СКФ такое разложение малоэффективно. Поэтому представляется более целесообразным построить разложение для другой линии, вдоль которой давление также равно давлению идеального газа (ЛИГ). Это линия для единичного фактора сжимаемости, представляющая собой в координатах плотность–температура прямую (для подавляющего числа веществ и модельных систем), проходящую в закритической области. Описание такого модифицированного вириального разложения приведено в [14]. В настоящей работе на основе этого подхода строится уравнение состояния для околокритической и сверхкритической областей. Это уравнение является единообразным для модельных систем и реальных веществ и определяется видом второго вириального коэффициента, критическими и бойлевскими параметрами.

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

ПРЕОБРАЗОВАНИЕ ВИРИАЛЬНОГО РЯДА

Рассмотрим систему частиц, взаимодействующих по известному парному потенциалу с характерными энергией ε и длиной r₀. Тогда плотность, температура, давление и s-й вириальный коэффициент (ВК) обезразмериваются с помощью соотношений n → $nr_{0}^{3}$, T → T/ε, P → ${{Pr_{0}^{3}} \mathord{\left/ {\vphantom {{Pr_{0}^{3}} \varepsilon }} \right. \kern-0em} \varepsilon }$, B_s → ${{{{B}_{s}}} \mathord{\left/ {\vphantom {{{{B}_{s}}} {r_{0}^{{3s - 3}}}}} \right. \kern-0em} {r_{0}^{{3s - 3}}}}$. Использование размерных единиц оговаривается отдельно.

Вириальное разложение, представляющее бесконечный ряд по степеням плотности, имеет вид [2]

Здесь P, n, Т, ${{B}_{s}}(T)$ – давление, плотность, температура и s-й вириальный коэффициент. Ряд (1) имеет плохую сходимость. Чтобы получить достаточно точное уравнение состояния, нужно учесть порядка десяти или даже более ВК. Для заданного потенциала взаимодействия ВК могут быть рассчитаны методом кластерных интегралов [2]. Однако этим методом относительно просто рассчитать второй и третий ВК. Коэффициенты более высокого порядка требуют существенно больших затрат компьютерного времени. Так, система частиц, взаимодействующих по потенциалу ЛД, используется для моделирования реальных веществ и для нее четвертый и пятый вириальные коэффициенты были вычислены более пятидесяти лет назад. Более старшие ВК появились совсем недавно [13]. В настоящее время вычислено 16 ВК для этой системы. Затем для коэффициентов, начиная с четвертого по девятый, была разработана аналитическая аппроксимация этих данных, справедливая в широком диапазоне температур [12, 13]. Для остальных ВК имеется только ряд температурных точек. Для других модельных систем ситуация во многом сходная. Методом кластерных интегралов рассчитано примерно 10 первых ВК, причем число доступных температурных точек уменьшается с увеличением номера коэффициента, а ошибка увеличивается. Результаты вычислений ВК для потенциала прямоугольной ямы представлены в [15, 16], для потенциала Морзе ‒ в [17], а для притягивательного потенциала Юкавы с твердым ядром ‒ в [18].

Результаты расчетов в вышеуказанных статьях показали, что ВК обнаруживают крайне нерегулярное поведение с изменением их номера. Это не позволяет установить какой-либо закономерности для прогнозирования их поведения при бóльших номерах s. Поэтому использование разложения (1), построенного при стремящейся к нулю плотности, малоэффективно для закритического флюида при увеличении плотности. Представляется более целесообразным использовать прямую линию для единичного фактора сжимаемости, лежащую в закритической области, вдоль которой давление равно давлению идеального газа [19]. Для многих веществ она обладает весьма интересным свойством, а именно для жидкостей, газов и сверхкритического состояния линия является абсолютно прямой на плоскости плотность–температура во всем диапазоне плотностей вплоть до линии кристаллизации. Более подробно о ее свойствах изложено в [10]. Здесь потребуется тот факт, что уравнение такой прямой линии можно записать в виде

где T_B и n_B – температура и плотность в точке пересечения с осью координат. Применение вириального разложения [10] показывает, что T_B ‒ это не что иное, как температура Бойля, т.е. B₂(T_B) = 0. Плотность n_B (далее бойлевская плотность) также выражается через второй и третий вириальные коэффициенты:

Значения n_B обычно близки к твердотельным [10], т.е. в этой области вириальный ряд уже может расходиться. Тем не менее ЛИГ продолжает выражаться уравнением (2) с бойлевскими параметрами, определяемыми вириальными коэффициентами.

Для дальнейшего удобно ввести эффективную температуру – аргумент ВК, в которой выделен комплекс (T/T_B + n/n_B), характеризующий ЛИГ. Эта температура определяется как

Для ЛД величины T_B = 3.418 и n_B = 1.14, а их отношение равно почти 3. Вдоль ЛИГ, когда справедливо уравнение (2), эффективная температура T* = T_B.

Для построения УРС с учетом ЛИГ используется тождество

Теперь правая часть уравнения (4) раскладывается в ряд по степеням х, а затем х полагается равным an. В результате для любого ВК получается ряд

Аналогичным образом раскладываются в ряды все вириальные коэффициенты, и получившиеся выражения подставляются в уравнение (1). Затем группируются слагаемые с одинаковыми степенями n. В результате получается обобщенный вириальный ряд вида

где модифицированные вириальные коэффициенты (МВК) D₂, D₃, D₄, … связаны с обычными B₂, B₃, B₄, … соотношениями

Ряд, состоящий из крайних правых членов в строках уравнений (6), содержащих B₂ и производные от B₂, суммируется

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

ВК для ЛД были получены в [13]. Для первых восьми из них построена надежная аналитическая аппроксимация, которая позволяет найти по уравнению (6) МВК для s = 2, 3, …, 8. Хотя МВК зависят и от плотности, и от температуры, что делает их сложней, чем ВК, их поведение оказывается более регулярным. Чтобы проиллюстрировать это, на рис. 1 представлена зависимость МВК как функция от B₂(T*).

Первое, на что следует обратить внимание, – МВК, в отличие от ВК, представляют гладкие функции. (Далее звездочка у температуры опускается там, где это не приводит к путанице.) Второе – при температуре Бойля, когда B₂(T_B) = 0, их значения становятся очень малыми. В масштабе рис. 1 при B₂(T_B) = 0 все кривые практически проходят через нулевую точку.

Коэффициент D₂(T_B) = B₂(T_B) = 0 по определению температуры Бойля. Равенство коэффициента D₃(T_B) = 0 вытекает из второй формулы уравнений (6) и уравнения (3). Действительно, после подстановки (3) в выражение для D₃(T_B) получается равенство

Значения ВК более высокого порядка в точке Бойля равны соответственно D₄ = –6.7 × 10^–4, D₅ = –1.6 × 10^–2, D₆ = 5.3 × 10^–2, D₇ = 5.3 × 10^–3, D₈ = 0.033.

При B₂(T) < B₂(T_B) с увеличением номера коэффициента соответствующие кривые сгущаются и стремятся к некоторой предельной кривой D_∞(T). Она была найдена в работе [14] и приведена на рис. 1 штриховой линией с символами.

При B₂(T) > B₂(T_B) МВК достигают максимума и затем начинают уменьшаться, стремясь к нулю. Для ВК с малыми номерами этот максимум лежит за рамками рис. 1. Однако с увеличением номера ВК максимум смещается в область малых температур, и он становится заметным, например, для D₇ и D_8.

По форме кривые на рис. 1 достаточно схожи с В₂. Более того, их можно совместить, смещая вверх и вправо на соответствующее расстояние. Это означает, что

где a₃, a₄, …, a_s и T₃, T₄, …, T_s – некоторые постоянные. Значения a₃, a₄, …, a_s находились из условия обращения в нуль МВК вдоль ЛИГ. Это дает

Далее величина смещения T_s находилась графически из условия наилучшего совпадения функций (6) и (7). Найденные таким образом значения T_s приведены в табл. 1.

Пунктиром на рис. 1 показаны коэффициенты D_s(T), рассчитанные по приближенным формулам (7). Как видно, аппроксимация (7) адекватно соответствует поведению первых восьми МВК. Однако ее можно экстраполировать на МВК более высокого порядка, если построить зависимость T_s = f(1/s). Такая зависимость представлена на рис. 2.

Зависимость T_s от s аппроксимируется прямой, уравнение которой имеет вид

Подставляя это уравнение в (7), получаем аппроксимацию МВК в виде

Отсюда легко найти значение МВК при s → ∞:

Данная зависимость показана на рис. 1 штриховой кривой с символами и соответствует линии D_∞(T), ранее полученной в работе [14].

Окончательно МВК принимает вид

Величина $\frac{A}{{{{T}_{{\text{B}}}}}} = \frac{{2.25}}{{3.418}} = 0.658$ в точности соответствует известному критическому инварианту [10] для ЛД, который равен

Как показано в работах [10, 11], для многих модельных систем и подавляющего числа реальных веществ величина I_c ≈ 0.6–0.8. Это дает основание переписать уравнение (8) в виде

где $I(T,n) = \frac{T}{{{{T}_{{\text{B}}}}}} + \frac{n}{{{{n}_{{\text{B}}}}}},$ ${{I}_{c}} = \frac{{{{T}_{c}}}}{{{{T}_{{\text{B}}}}}} + \frac{{{{n}_{c}}}}{{{{n}_{{\text{B}}}}}}$.

Из (10) следует, что вдоль ЛИГ, когда I(T, n) = 1, давление равно давлению идеального газа.

Уравнение (10) справедливо, если аргументы у В₂ положительны, т.е.

На плоскости n–T это неравенство описывает семейство прямых, параллельных линии I(T, n) = 1 (ЛИГ) и смещенных в область меньших значений температуры и плотности. Величина смещения меняется от нуля для s = 2 до 0.658 при s → ∞. Тем самым возникает ограничение применимости настоящего подхода для малых плотностей и температур. Однако это ограничение возникает вне области критических и сверхкритических флюидов и поэтому не представляет интереса для данной работы.

После всех преобразований давление имеет вид

Вдоль ЛИГ I(T, n) = 1, поэтому каждое слагаемое в сумме обращается в нуль. Давление становится равным давлению идеального газа. В критической точке МВК принимают вид

Учитывая (9), правую часть (13) можно переписать как ${{B}_{2}}\left( {{{2{{T}_{{\text{B}}}}0.658} \mathord{\left/ {\vphantom {{2{{T}_{{\text{B}}}}0.658} s}} \right. \kern-0em} s}} \right) - {{B}_{2}}({{2{{T}_{{\text{B}}}}{\text{(0}}{\text{.171}}\,\,{\text{ + }}\,\,{\text{0}}{\text{.658}}} \mathord{\left/ {\vphantom {{2{{T}_{{\text{B}}}}{\text{(0}}{\text{.171}}\,\,{\text{ + }}\,\,{\text{0}}{\text{.658}}} s}} \right. \kern-0em} s}))$.

Нужно заметить, что в (12) из дополнительных параметров входит только функциональная форма B₂(T). Она же определяет значение T_B. Сам же коэффициент B₂(T) легко вычисляется для любой системы с заданным парным потенциалом взаимодействия как интеграл от функции Майера [2]. Поэтому сначала следует проверить (12) именно для таких систем.

СРАВНЕНИЕ С РЕЗУЛЬТАТАМИ РАСЧЕТОВ ДЛЯ СИСТЕМЫ ЛЕННАРДА–ДЖОНСА

Система ЛД является основной тестовой системой для численных расчетов в статистической физике классических систем уже более полувека [2]. Поэтому для ЛД также построены референтные уравнения состояния аналогично тому, как это было сделано для реальных веществ [4]. Одно из наиболее точных и недавних таких РУРС, полностью описывающее численные данные по термодинамике ЛД, было построено в [20]. Оно использовалось в настоящей работе для проверки уравнения (12). Критические параметры для ЛД указаны в табл. 2. На рис. 3 представлены изотермы давления, полученные в данной работе согласно уравнению (12), и изотермы из работы [20] в приведенных к критическим координатах.

Как видно, сверхкритические изотермы с T = 1.5 и 2.5 хорошо согласуются с данными [20]. Критическая изотерма данной работы несколько расходится с численными данными в области высоких плотностей. При низких плотностях она обрывается в силу нарушения неравенства (11).

СРАВНЕНИЕ С ДРУГИМИ МОДЕЛЬНЫМИ СИСТЕМАМИ

Рассмотрим теперь как уравнение (12) работает для систем с потенциалами, отличными от ЛД. В качестве таковых рассмотрим упомянутые выше потенциалы прямоугольной ямы и Юкавы с твердым ядром.

Потенциал прямоугольной ямы имеет вид

где σ – диаметр твердой сердцевины; λ – параметр, характеризующий протяженность потенциальной ямы; ε – глубина ямы притяжения. Второй и третий вириальные коэффициенты вычисляются аналитически и приведены, например, в [21]. Все они зависят от параметра λ. В [22] было отмечено, что при слишком малых или слишком больших λ линия единичного фактора сжимаемости искривляется. Здесь выбрано значение λ = 1.75, при котором эта линия является прямой.

Рассматриваемый потенциал Юкавы имеет вид

При k = 1.8 ЛИГ остается прямой, как показали расчеты в [22]. Там же были рассчитаны и приведены бойлевские параметры для целого ряда модельных потенциалов, включая (14), (15). Координаты критической точки для потенциала прямоугольной ямы были численно рассчитаны в [23], а для потенциала Юкавы – в [24]. Эти параметры представлены в табл. 2.

Результаты расчета критических изотерм и критических точек приведены на рис. 4 и в табл. 2. Как видно, положение критических точек по результатам численного моделирования находится в хорошем соответствии с рассчитанными в данной работе критическими изотермами. Кроме того, в табл. 2 представлены значения ΔH – безразмерная энтальпия парообразования в тройной точке, деленная на число частиц, Z_c – фактор сжимаемости в критической точке, а также ряд безразмерных параметров, характеризующих асимметрию бинодали и используемых в различных законах подобия, (см., например, [25]). Для ЛД ΔH рассчитывается прямо по РУРС [20]. Для других двух систем ΔH не представлена в явном виде, но в [23, 24] есть зависимости давления P и разности удельных объемов фаз Δν от температуры Т на бинодали. Это позволило получить ΔH по уравнению Клапейрона–Клаузиуса [2], как ΔH = TP 'Δν, где P ' = dP/dT вдоль бинодали. Значения температур в тройной точке Т_tr для этих систем: Т_tr = 0.661 – ЛД [20], Т_tr = 0.9 – потенциал Юкавы [24], Т_tr = 1.0 – прямоугольнaя ямa [23].

СРАВНЕНИЕ С КРИТИЧЕСКИМИ ПАРАМЕТРАМИ МЕТАНА

Рассмотрим теперь, как работает уравнение (12) для реальной системы. В отличие от модельных систем здесь нет явно заданного потенциала взаимодействия, однако для многих веществ второй вириальный коэффициент может быть определен из экспериментальных данных при низких плотностях [2]. Благородные газы в рассматриваемой области неплохо описываются потенциалом ЛД. И для них можно ожидать, что уравнение (1) будет вести себя, как и для ЛД. Поэтому интересно рассмотреть вещество, не описываемое потенциалом ЛД, но для которого есть достаточно точные данные измерений и ЛИГ – прямая. Таковым был выбран метан CH₄. Для метана B₂(T) в температурном диапазоне 115–670 К приведен в [26], где для него построена аппроксимирующая формула. Эти данные позволяют найти температуру Бойля Т_B = 509.4, что практически совпадает со значением 510 К, которое было получено экстраполяцией ЛИГ к нулевой температуре [11]. Значения критических параметров для CH₄ таковы Т_с = 190.6 К, n_c = 0.163 г/cм³, P_c = 46 атм [4]. Полагая в (9) I_c = 0.67, как и для многих других веществ, получаем n_B = 0.55 г/см³, что тоже очень мало отличается от значения 0.57 г/см³, полученного экстраполяцией ЛИГ к нулевым плотностям [11]. Уравнение (12) в абсолютных единицах принимает вид

где k – постоянная Больцмана, М – масса атома (для метана M = 6.64 × 10^–24 г).

Фактор сжимаемости

В (17) введена плотность b, связанная с исключенным объемом атома. В ЛД эта область определяется расстоянием, на котором отталкивательная ветвь потенциала пересекает ось абсцисс. Для реальных веществ она зависит от размера атома и по порядку величины равна b ~ М/R³, где R – размер атома. Ниже используется значение b = 0.87 г/см³, чему соответствует R ~ 2 Å. Эта величина несколько больше длины углерод-водородной связи в молекуле метана, равной 1.54 Å.

На рис. 5 по формуле (16) построена критическая изотерма 2 для метана c вышеперечисленными параметрами и значением b = 0.87 г/см³. Соответствующие критические параметры, полученные в подходе данной работы: Т_с = 208.9 К, n_c = = 0.143 г/cм³, P_c = 40 атм. Здесь же нанесена соответствующая экспериментальная изотерма 1, построенная по данным [4]. Символами отмечены положения критических точек. Как видно, имеет место согласие между кривыми 1 и 2. Расхождение лежит, по-видимому, в пределах точности методов, которые использовались при расчете этих кривых.

ЗАКЛЮЧЕНИЕ

Построено единообразное уравнение состояния в околокритической или в закритической области. Оно было использовано для реального вещества (метана) и трех модельных систем. Это уравнение включает в себя физически обоснованные и определяемые видом потенциала величины (второй вириальный коэффициент, бойлевские и критические параметры) и не содержит каких-либо других эмпирических констант. Формально применимость рассматриваемого уравнения ограничена со стороны области c низкой плотностью, неудовлетворяющей условию (11). Однако эта область по свойствам близка к идеальному газу и не представляет особого интереса для критических или сверхкритических флюидов. Критические параметры, рассчитанные с помощью рассматриваемого уравнения, находятся в согласии с литературными данными. Заметим также, что интегрирование уравнения (12) по плотности позволяет получить выражение для свободной энергии, а дифференцирование последней по температуре – выражения для внутренней энергии и энтропии [2]. Обычные вириальные коэффициенты зависят только от температуры, поэтому для обычного вириального ряда такое интегрирование осуществляется аналитически. Но МВК зависят как от температуры, так и от плотности, поэтому для (12) данная процедура не приведет к простому аналитическому выражению. Тем не менее подобное исследование можно провести в дальнейшем. Кроме того, планируется рассмотреть, как предложенное уравнение ведет себя и для других модельных систем. Это могут быть те же системы, но в двухмерном пространстве, для которых линия Z = 1 тоже прямая [27], или модифицированная однокомпонентная модель плазмы классических ионов на однородно-сжимаемом фоне идеального газа электронов, рассмотренная в [28].

Работа выполнена при финансовой поддержке Министерства науки и высшего образования РФ (соглашение с ОИВТ РАН № 075-15-2020-785 от 23 сентября 2020 г.).