Теплофизика высоких температур, 2022, T. 60, № 5, стр. 682-691

ВВЕДЕНИЕ

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

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

Теоретические модели, основанные на потенциалах взаимодействия молекул [6–29], позволяют рассчитывать термодинамические свойства плотного газа в согласии с экспериментальными данными и результатами моделирования методами Монте-Карло и молекулярной динамики. Совершенствование аппарата статистической механики стимулировало разработку вариационной теории MCRSR (Mansoori–Canfield–Rasaiah–Stell–Ross) [9, 10], термодинамических теорий возмущений [6, 7, 11–13] и интегральных уравнений для функций распределения молекул в приближении мягких сфер HMSA (hypernetted-chain/soft core mean spherical approximation) [14, 28].

Использование теории возмущений [6, 7, 11–13, 15, 25–27, 29] для разработки УРС плотных газов основано на предположении, что структура флюида в первую очередь определяется силами отталкивания на малых расстояниях и лишь в незначительной степени силами притяжения большого радиуса действия. Поэтому в теории возмущений для сферических молекул в качестве базисной используется система твердых сфер, а слабые силы притяжения учитываются в виде возмущений.

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

В настоящей работе каждая методика разделения потенциала расширена для случая бинарных смесей. Проведена верификация результатов термодинамического моделирования свойств бинарных смесей (H₂–He, NH₃–H₂, N₂–N, O₂–O) с доступными экспериментальными данными и результатами компьютерного моделирования методами Монте-Карло и молекулярной динамики. Для проведения термодинамических расчетов использовалась разработанная авторами теоретическая модель уравнения состояния двухкомпонентной системы [15]. Впервые выполняется сравнение предсказаний теоретической модели УРС двухкомпонентной смеси [15] с использованием различных версий теории возмущений с результатами моделирования Монте-Карло и экспериментальными данными термодинамических свойств бинарных смесей.

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

ОСНОВНЫЕ ПОЛОЖЕНИЯ ТЕОРИИ ВОЗМУЩЕНИЙ

Основная идея любой теории возмущений состоит в разделении потенциала φ_ij(r) исследуемой системы в точке r = λ на две составляющие: основную часть $\varphi _{{ij}}^{{{\text{ref}}}}(r),$ представляющую собой потенциал взаимодействия молекул в некоторой базисной системе (ref – reference – базисный), и малое возмущение $\varphi _{{ij}}^{{{\text{pert}}}}(r)$ (pert – perturbation – возмущение). В качестве базисной системы обычно используется модель твердых сфер [10–15], так как она наиболее изучена и для нее известна радиальная функция распределения молекул.

Для двухкомпонентной смеси разделение потенциала взаимодействия φ_ij(r) для каждой пары молекул ij (ij = 1, 2) можно представить в виде

где r – расстояние между молекулами, λ_ij – точка разделения потенциала, φ_λ,ij(r) – сглаживающий потенциал.

Методика расчета базисной и возмущающей частей потенциала (1), а также определение точки разделения потенциала целиком определяются конкретной версией теории возмущений.

Для расчета термодинамических свойств равновесной системы при заданных температуре T и объеме V достаточно знать энергию Гельмгольца этой системы. Тогда величины всех необходимых термодинамических параметров могут быть найдены путем вычисления соответствующих производных энергии Гельмгольца.

Для двухкомпонентной системы избыточная свободная энергия Гельмгольца F^ex в расчете на одну молекулу в соответствии с теорией возмущений определяется по формуле

где x_i, x_j – мольные доли компонентов смеси; ρ – плотность смеси; $g_{{ij}}^{{HS}}\left( {r,\,\,{{\rho }},{\text{\;}}T} \right)$ – радиальная функция распределения твердых сфер, зависящая от расстояния между молекулами типа i и j; β = 1/(k_BT); k_B – постоянная Больцмана, F ^HS – избыточная энергия Гельмгольца базисной системы твердых сфер, методика расчета которой описана в [16].

Для расчета радиальных функций распределения $g_{{ij}}^{{HS}}\left( {r,\,\,{{\rho }},{\text{\;}}T} \right)$ в [17] применялась методика, включающая корректирующую поправку, предложенную Верле и Вейсом [18]:

где $g_{{ij}}^{{{\text{PY}}}}$(r) – функции распределения твердых сфер, полученные из уравнения Перкуса–Йевика [18]; R_ii – эффективный диаметр твердой сферы i-го компонента; $\xi = \frac{\pi }{6}\sum\nolimits_{i = 1}^2 {{{\rho }_{i}}R_{{ii}}^{3}} $ – коэффициент упаковки твердых сфер; $R_{{ij}}^{w}$ определяется по формуле $R_{{ij}}^{w}$ = R_ij(1 – ξ/16)^1/3.

Использование выражения (3) позволяет получить более точные значения радиальных функций распределения $g_{{ij}}^{{HS}}\left( {r,{{\rho }},{\text{\;}}T} \right)$, чем функций $g_{{ij}}^{{{\text{PY}}}}$(r), полученных из уравнения Перкуса–Йевика [19], в сравнении с данными моделирования методом Монте-Карло. Аналитическое выражение для расчета функций распределения, а также подробное описание методики представлено в [17].

Численный расчет интеграла, присутствующего в формуле (2), требует значительных временных затрат, быстро увеличивающихся с ростом переменной интегрирования r. Поэтому в настоящей работе используется методика, которая без потери точности значительно упрощает вычисления и сокращает время расчета. Второе слагаемое формулы (2), содержащее интеграл для i, j = 1, 2, с учетом выражения (3) и разделения потенциала (1) может быть представлено в виде суммы трех составляющих:

Появление слагаемого F_3,ij в выражении для энергии Гельмгольца (4) связано с коррекцией (3) функций распределения твердых сфер Перкуса–Йевика согласно методике Верле и Вейса [18]. В [29] показано, что эта коррекция значительна только при высоких плотностях. Поэтому вклад в суммарную энергию Гельмгольца оказывается наибольшим именно в области высоких плотностей. Авторы [27] в своей методике расчетов по теории KLRR пренебрегают слагаемым (7) в формуле (4), безосновательно указывая, что вклад в суммарную величину энергии Гельмгольца невелик.

Также для численного вычисления интегралов (5)–(7) использовалась методика, подробно описанная в [15] и позволяющая рассчитать избыточную энергию Гельмгольца с высокой точностью и быстродействием по сравнению с прямым численным интегрированием в выражении (3).

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

Например, коэффициент сжимаемости Z = =  PV/(Nk_BT), соответствующий давлению P, рассчитывается путем численного дифференцирования избыточной энергии Гельмгольца по плотности:

Аналогично вычисляется избыточная внутренняя энергия в расчете на одну молекулу путем дифференцирования по температуре:

Входящие в выражения (8) и (9) частные производные избыточной энергии Гельмгольца (4) по переменным ρ и β вычисляются численным дифференцированием функции на равномерной сетке по четырем соседним узлам. Для этого применяется формула, полученная по методу Рунге–Ромберга, обеспечивающая достаточный для целей данной работы четвертый порядок точности по шагу сетки. Эта методика позволяет избежать появления нефизичных слагаемых при дифференцировании избыточной энергии Гельмгольца при использовании термодинамически зависимых потенциалов взаимодействия. Через соответствующие первые и вторые производные энергии Гельмгольца также можно выразить и все остальные термодинамические параметры флюида.

Теория возмущений Баркера–Хендерсона. Для случая бинарных смесей в версии Баркера и Хендерсона (BH) [12] базисный потенциал φ^ref(r) и возмущающая часть φ^pert(r) имеют вид

Здесь r – межмолекулярное расстояние; σ_ij – расстояние, при котором φ_ij(σ_ij) = 0, является точкой разделения межмолекулярного потенциала взаимодействия φ_ij(r).

Версия BH успешно использовалась для расчета термодинамических параметров однокомпонентных систем, молекулы которых взаимодействуют в соответствии с потенциалом Леннард–Джонса [20], а также при расчете равновесных параметров системы пар–жидкость при высоких давлениях [21].

Теория возмущений Викса, Чендлера, Андерсена (WCA). Теория WCA [13, 22] – широко используемая модель теории возмущений, в которой базисный потенциал φ^ref(r) и возмущающая часть φ^pert(r) определяются как

Отличие WCA от BH заключается в определении базисного потенциала и местоположения точки разделения потенциала взаимодействия. В версии WCA точка разделения потенциалов совпадает с местоположением минимума потенциала межмолекулярного взаимодействия r_m,ij = 2^1/6σ_ij.

Для определения диаметров твердых сфер предложен термодинамический критерий – равенство сжимаемостей базисной системы и эквивалентного твердосферного флюида. Эта методика позволяет описать термодинамические свойства плотных однокомпонентных флюидов в хорошем согласии с экспериментальными данными [13, 22, 23] лишь в области низких плотностей, хотя результаты могут быть реалистичными и при высоких плотностях в зависимости от вида потенциала межмолекулярного взаимодействия.

На основе метода WCA в [22] получено выражение для энтропии, которое успешно используется для расчета термодинамических параметров различных сплавов.

Однако результаты расчетов термодинамических свойств однокомпонентных веществ с использованием версий теории возмущений BH и WCA отклоняются от результатов моделирования методом Монте-Карло на 10–15% и дают нереалистичные результаты в области высоких плотностей и температур, представляющих наибольший интерес.

Версия KLRR теории возмущений. Более удачной из теорий, позволяющих получать УРС плотных газов как при высоких давлениях и температурах, так и при более низких температурах и высоких плотностях, является теория возмущений KLRR [11]. Эта теория воспроизводит результаты расчетов Монте-Карло с хорошей точностью для разных типов потенциалов взаимодействия молекул, включая потенциалы exp-6, реалистичность которых для расчетов параметров ударного сжатия и детонации подтверждается многими исследованиями [25, 26].

В версии KLRR сглаживающий потенциал φ_λ(r) является функцией межмолекулярного расстояния, в отличие от предыдущих версий (нулевое значение в BH и постоянная величина в WCA):

Определение точки разделения потенциала λ_ij и функциональной формы сглаживающего потенциала φ_λ,ij(r) является произвольным при выполнении условия, что $\varphi _{{ij}}^{{{\text{ref}}}}$(r) имеет отталкивательный характер.

В данной работе используется самая простая – линейная – форма для сглаживающего потенциала φ_λ,ij(r):

Точка разделения потенциала межмолекулярного взаимодействия в оригинальной версии KLRR зависит от плотности и вычисляется как

где a_fcc,i = 2^1/6/$\rho _{i}^{{1{\text{/}}3}}$ – эквивалент расстояния ГЦК-решетки, ρ_i = N_i/V – концентрация i-го компонента в смеси.

Эта схема переходит в форму WCA при низких плотностях, имеющих параметр решетки a_fcc,i > r_m,ii. Для плотностей выше кривой плавления и температур выше k_BT/ε_ii = 100 использование этой схемы приводит к согласию результатов расчетов с данными метода Монте-Карло для однокомпонентных флюидов, молекулы которых взаимодействуют по потенциалам Леннард–Джонса, exp-6 и обратностепенным потенциалом.

Однако метод расчета точки разделения λ_ii имеет недостаток, связанный с выбором негладкой функции (13). Скачкообразная зависимость функции λ_ii(ρ_i) приводит к тому, что вторые и высшие частные производные избыточной энергии Гельмгольца по плотности и термодинамические параметры, вычисляемые по этим производным (например, коэффициент изотермического сжатия), имеют разрывы при той же плотности. Этот недостаток теории KLRR может приводить к практическим трудностям при выполнении термодинамических расчетов с таким УРС.

В настоящей работе значения точек разделения для одноименных потенциалов λ_ii вычисляются с использованием следующей непрерывной функции, предложенной в [27] для однокомпонентного флюида:

где ρ_σ,ii = ρ_ii/√2, n – некоторый параметр.

При n = 120 функции (13), (14) становятся очень близки друг к другу. На рис. 1 показана зависимость точки разделения λ_ii от плотности при различных значениях параметра n в формуле (14).

Таким образом, замена кусочно заданной функции (13) в правиле выбора точки разбиения потенциала на гладкую функцию (14) с n = 120 позволяет устранить недостаток оригинальной версии KLRR [11] и тем самым избежать разрывов в производных эффективного диаметра твердых сфер.

Модифицированная таким образом версия теории возмущений KLRR [11] с учетом выражений для расчета энергии Гельгольца (5)–(7) является улучшенной версией – KLRR-T, на основе которой разработана теоретическая модель УРС [15], реализованная в виде вычислительных программ, для расчета термодинамических свойств бинарных смесей. При этом теория сохраняет свои главные достоинства: высокую точность в широком диапазоне изменений плотности и температуры и быстродействие.

С использованием (14) в данной работе рассчитываются местоположения разделения потенциалов лишь одноименных взаимодействий (i = j = 1, 2 для бинарной смеси). Точка разделения λ_ij перекрестного потенциала взаимодействия определяется таким образом, чтобы выполнялось условие аддитивности диаметров твердых сфер базисной системы R₁₂ = (R₁₁ + R₂₂)/2.

РАСЧЕТ ТЕРМОДИНАМИЧЕСКИХ ПАРАМЕТРОВ БИНАРНЫХ СМЕСЕЙ С ПОТЕНЦИАЛОМ ВЗАИМОДЕЙСТВИЯ EXP-6

Для тестирования рассмотренных версий теории возмущений для расчета термодинамических параметров бинарных смесей взаимодействие между частицами в системе описывается модифицированным потенциалом Букингема exp-6

где ε_ij > 0 – глубина потенциальной ямы (min[φ_ij(r)] = = −ε_ij); r_m,ij – расстояние между молекулами, при котором потенциальная энергия минимальна (φ_ij(r_m,ij) = −ε_ij); α_ij – параметр, определяющий жесткость отталкивания молекул.

Однако трехпараметрический потенциал (15), (16) с постоянными значениями параметров α_ij, r_m,ij, ε_ij хорошо описывает межмолекулярные силы только для веществ с неполярными или почти неполярными молекулами, но оказывается слишком грубым приближением для полярных молекул, в особенности для молекул с водородными связями (например, H₂O и NH₃) [30].

В [30] на примере H₂O показано, что зависящие от ориентации электростатические взаимодействия молекул могут быть учтены в квантово-механических потенциалах путем введения в эти потенциалы температурных зависимостей. В [27] приведена соответствующая модифицированная форма и для потенциала exp-6, в которой глубина потенциальной ямы зависит от температуры:

где θ_ii – параметр, отвечающий за учет электростатических эффектов.

Дополнение потенциала температурной зависимостью позволяет учитывать вклад от дипольных взаимодействий, который становится значительным при температурах порядка θ_ii и ниже. Модифицированный таким образом потенциал exp-6 (15)–(17) имеет четыре параметра α_ii, r_m,ii, ε_0,ii, θ_ii, причем для неполярных веществ θ_ii = 0. На примере воды и аммиака показано [25, 26, 30, 31], что термодинамический расчет с теоретическим УРС на основе потенциалов вида (15)–(17) удовлетворительно воспроизводит результаты ударно-волновых измерений удельного объема этих веществ в широком диапазоне изменения давления (от тысяч до сотен тысяч атмосфер). При использовании стандартного трехпараметрического потенциала exp-6 (15)–(17) такие результаты получить не удается.

Методика определения параметров потенциалов exp-6 (15)–(17) по данным динамических (ударно-волновых) и статических экспериментов разработана в [29]. В настоящей работе для дальнейших расчетов используются потенциальные параметры для одноименных пар молекул, найденные в [29], которые представлены в табл. 1.

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

Здесь k_ij, l_ij и m_ij – корректирующие множители, значения которых обычно близки к единице. Эти множители могут быть определены, если имеется подходящая теоретическая информация или экспериментальные данные. В противном случае принимается k_ij = l_ij = m_ij = 1, что сводит формулы (18) к классическим смесевым правилам Лоренца–Бертло. В данной работе для описания перекрестных взаимодействий используются аддитивные параметры, т.е. k_ij = l_ij = m_ij = 1 в выражениях (18).

Таким образом, в предложенной модели УРС двухкомпонентной смеси плотных газов межмолекулярные взаимодействия описываются модифицированными, зависящими от температуры потенциалами exp-6 (15)–(17). Тем самым разработанная модель УРС оказывается применимой для смесей, содержащих вещества как с неполярными, так и с полярными молекулами.

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

Термодинамические расчеты выполнены с использованием модели уравнения состояния бинарных смесей [15]. Модель УРС [15] реализована в виде вычислительного алгоритма и компьютерной программы, позволяющей рассчитывать избыточные термодинамические свойства флюида при заданных значениях температуры и плотности. Все исследуемые в работе версии теории возмущений реализованы в этой программе путем редактирования соответствующих блоков подпрограмм, в которых производится разделение потенциала взаимодействия на базисную и возмущающую части в соответствии с уравнениями (10), (11) или (12).

На рис. 2 показан вид возмущающей части потенциала и местоположения точек разделения потенциала для каждой версии теории возмущений на примере молекулы азота.

Для тестирования первой выбрана бинарная смесь гелия и водорода H₂–He. Выбор этой смеси объясняется наличием данных моделирования методом Монте-Карло [32] в широком интервале температур 50–7000 К и плотностей ρ = 0.15–0.8 г/см³, а также значительно различающимися потенциальными параметрами компонентов смеси. Однако в [32] расчеты выполнены с использованием неаддитивного перекрестного потенциала для H₂–He. Поэтому дополнительно проведено моделирование методом Монте-Карло с помощью программного кода MC Towhee [33, 34] в указанном интервале термодинамических состояний.

Результаты расчетов представлены в табл. 2, в которой показаны исходные параметры (температуры, плотность, состав смеси), рассчитанные значения давления и внутренней энергии с использованием различных схем разделения потенциала взаимодействия, а также результаты моделирования методом Монте-Карло (МК).

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

Наилучшее согласие с данными моделирования методом Монте-Карло имеют результаты расчетов, полученные с использованием версии KLRR-T теории возмущений, особенно в области высоких температур и давлений (табл. 2, 3).

Далее рассматриваются бинарные смеси N₂–N и O₂–O, которые являются продуктами ударно-волнового сжатия жидких азота и кислорода вследствие диссоциации при высоких температурах. Ударные адиабаты N₂ и O₂ тщательно исследованы, доступно большое количество экспериментальных и теоретических результатов для сравнения [35‒42].

Моделирование термодинамических параметров бинарных смесей N₂–N и O₂–O, полученных в результате ударно-волнового сжатия, также выполнено с использованием всех трех способов разделения межмолекулярного потенциала взаимодействия. Результаты расчетов в соответствии с экспериментальными данными и результатами моделирования МК [35, 44] представлены графически в виде ударных адиабат в P–V- и T–P-координатах на рис. 3, 4. Следует отметить, что на рис. 3б заметно отклонение результатов расчета по KLRR-T от экспериментальных данных [35] при давлении выше 50 ГПа. Это связано с началом процесса ионизации в этой области, которая не учитывалась в авторских расчетах для исследуемых версий теории возмущений. Поэтому дополнительно выполнено сравнение с данными МК-моделирования [35, 44] во всем рассматриваемом диапазоне термодинамических состояний. В табл. 4 представлена статистика относительных отклонений значений давления на ударной адиабате жидких азота и кислорода от данных МК-моделирования [43] для количественной оценки точности рассчитываемых свойств бинарных смесей N₂–N и O₂–O на основе исследуемых в работе версий теории возмущений.

Рис. 3, 4 и табл. 4 демонстрируют, что результаты термодинамического моделирования с использованием версии KLRR-T [15] теории возмущений для расчета состояний ударно-волнового сжатия жидких азота и кислорода согласуются с экспериментальными данными в пределах указанной авторами погрешности, расчетами других авторов и результатами моделирования методом Монте-Карло. В свою очередь версии BH [12] и WCA [13] лишь удовлетворительно описывают свойства бинарных смесей N₂–N и O₂–O, образующихся в результате диссоциации при сжатии в ударной волне.

В [45] опубликованы экспериментальные данные для бинарных смесей, одним из компонентов которых является аммиак. Эти данные предоставляют хорошую возможность верифицировать описанные методики разделения потенциалов для случая двухкомпонентных смесей, содержащих полярные молекулы аммиака. Это смеси аммиак–водород, содержащие 34.75 и 57.49% аммиака. Результаты расчетов показаны графически в виде изотерм при 423 и 573 К в P–V-координатах на рис. 5. Статистика отклонений результатов расчетов от экспериментальных данных представлена в табл. 5.

Рис. 5 и табл. 5 демонстрируют одинаково хорошее согласие с экспериментальными данными результатов на основе версий KLRR-T и WCA.

Анализируя полученные результаты, можно отметить, что погрешность расчетов на основе версий BH [12] и WCA [13] растет с увеличением отношений потенциальных параметров компонентов бинарных смесей. Это поведение объясняется тем, что способы разделения потенциалов взаимодействия и точки разделения потенциалов напрямую связаны со значениями параметра r_m и не зависят от полярности молекул веществ смеси. Таким образом, можно предположить, что эти версии теории возмущений могут применяться только для моделирования свойств бинарных систем с незначительно различающимися потенциальными параметрами r_m.

ЗАКЛЮЧЕНИЕ

В результате рассмотрения различных версий теории возмущений показано, что версия теории возмущений KLRR-T предсказывает термодинамические свойства бинарных смесей веществ с заметно различающимися потенциальными параметрами в согласии с экспериментальными данными и результатами моделирования методами Монте-Карло и молекулярной динамики в широкой области высоких давлений и температур. Эта область охватывает как низкоплотные состояния разреженного газа, так и высокоплотные состояния с давлениями в десятки гигапаскалей, пока влияние ионизации не начнет изменять параметры состояния продуктов ударно-волнового сжатия веществ.

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