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

ВВЕДЕНИЕ

Задача построения аналитических уравнений состояния (УС), позволяющих с удовлетворительной точностью прогнозировать свойства вещества в широком диапазоне изменения параметров состояния, до сих пор привлекает к себе постоянное внимание.

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

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

Теоретическое исследование структуры и свойств жидкостей может быть проведено на основе континуальной модели, состоящей из N одинаковых частиц, заключенных в объеме V и взаимодействующих согласно потенциалу v(r).

Исследования этой модели, выполненные в ряде работ [4–6] с привлечением различных комбинаций потенциалов Юкавы (в том числе двойной – DY-потенциал), потенциала Леннард-Джонса и др., где применяется техника интегральных уравнений и методы численного моделирования, содержат результаты, позволяющие провести сравнение с экспериментом. Рассмотрение в [7, 8] трехмерной модели с обобщенным потенциалом Морзе с помощью уравнения Перкуса–Йевика ограничивается построением уравнения состояния, для которого не приводятся конкретные расчетные данные, что затрудняет оценку пригодности УС для описания простых жидкостей. Здесь, как и в [7, 8], под обобщенным потенциалом Морзе следует понимать линейную комбинацию двух экспоненциальных потенциалов вида C exp(–cr). В специальной литературе, например в [9], можно встретить и другое название — двойной экспоненциальный потенциал (далее DE-потенциал).

Вопрос о предсказательных возможностях модели с DE-потенциалом, по-видимому, не потерял своей актуальности и может быть частично разрешен в рамках подхода, предложенного ранее в [10], где для системы с парным изотропным потенциалом v(|r|) было получено следующее выражение для свободной энергии Гельмгольца:

Здесь β = 1/(k_BT) – обратная температура, k_B – постоянная Больцмана, n = N/V – концентрация частиц, F_id= N(ln(nλ³) − 1)/β, $\lambda = {h \mathord{\left/ {\vphantom {h {\sqrt {2\pi {{m}_{0}}{{k}_{{\text{B}}}}T} }}} \right. \kern-0em} {\sqrt {2\pi {{m}_{0}}{{k}_{{\text{B}}}}T} }}$ – тепловая длина волны де Бройля, h – постоянная Планка. Предполагается, что потенциал v(|r|) имеет фурье-образ ${\tilde {v}}(k).$

Важно подчеркнуть, что в [10] формула (1) получена путем сведения КИ к интегралу типа Лапласа, вычисленного в квадратичном приближении метода перевала, а в пионерской работе Д.Н. Зубарева [11] аналогичная формула выведена с привлечением метода коллективных переменных. Кроме того, этот результат следует из выражения для КИ, найденного в [12] с помощью эргодической теоремы Вейля.

Перспективы использования выражения (1) для расчета свойств модельных систем открываются после вычисления интеграла I(n, β). В случае простейших двухпараметрических потенциалов [13, 14] интеграл (2) легко вычисляется. Однако практический интерес представляют “реалистичные” потенциалы, содержащие более двух регулируемых параметров, варьированием которых можно попытаться улучшить согласие теоретических результатов с экспериментом. По мере усложнения потенциала задача нахождения I(n, β) становится более трудоемкой, а выражение (1) менее приспособленным для проведения расчетов.

Преодоление возникающих на этом пути трудностей может быть связано, по крайней мере, с двумя подходами:

– если найдено точное, но достаточно громоздкое выражение для I(n, β), то его можно адаптировать для удобных аналитических расчетов, привлекая соображения эвристического характера или рассматривая предельные случаи [15, 16];

– если же точное выражение для I(n, β) найти невозможно либо оно оказывается необозримым и применение выражения (1) для дальнейших расчетов теряет смысл, то для вычисления (2) необходимо применять теоретически обоснованные приближенные методы.

Следуя первому из указанных подходов, удается получить адекватные результаты, относящиеся к системам с четырехпараметрическим осциллирующим потенциалом [15] и двойным потенциалом Юкавы (DY-потенциал) [16]

который выступает здесь в роли “пробного камня” с целью проверки эффективности второго подхода, используемого для получения УС системы с DE-потенциалом: Параметры a, b, A, B положительны. Применение потенциалов для расчетов свойств в рамках приближения (1) связано с неравенством ${\tilde {v}}(k)$ ≥ 0 [17, 18], определяющим область устойчивости взаимодействий, где невозможен коллапс системы. Решение этого неравенства с фурье-образами потенциалов (3) и (4) (w_DY= A/a², t = k/a, w_DE= A/a⁴) в системе координат δ = b/a и ε = B/A представляет собой фигуру, образованную полосой G₂= {(δ, ε): δ ≥ 1, 0 ≤ ε ≤ 1} и сектором G₁= {(δ, ε): 0 < δ < 1, 0 < ε < δ^m} с m = 2 для потенциала (3), m = 4 для (4). Области устойчивости изображены на рис. 1. Значения δ и ε в полосе G₂ таковы, что при любых r в (3) и (4) лидирует первое слагаемое и потенциалы преобразуются к чисто отталкивательным. Поэтому интерес представляет сектор G₁, где параметры δ и ε принимают такие значения, при которых потенциалы (3), (4) имеют черты, характерные для “реальных” непрерывных потенциальных функций: отталкивание, притяжение и минимум, как следствие, первых двух свойств.

УРАВНЕНИЯ СОСТОЯНИЯ

При интегрировании логарифма ln(1 + nβ${\tilde {v}}(t)$) в (2) с фурье-образом (6) возникает дробно-рациональная функция со знаменателем в виде многочлена высокой степени относительно t, что делает весьма затруднительным дальнейший поиск точного выражения для I(n, β), которое будет иметь настолько громоздкий вид, что вряд ли может быть размещено на одной странице журнала, даже с учетом вспомогательных обозначений.

Для устранения таких технических сложностей предлагается простой способ преобразования подынтегрального логарифма. Реализация этого способа сводится к такому представлению аргумента 1 + nβ${\tilde {v}}(t)$ в виде произведения двух сомножителей, при котором соответствующая сумма логарифмов содержит слагаемое, определяющее основной вклад в интеграл и интегрируемое без каких-либо приближений. Тогда, пренебрегая малым слагаемым и сохраняя под знаком интеграла доминирующее, можно получить приближенное выражение для I(n, β). Для проверки работоспособности предлагаемой вычислительной схемы привлекается модель с DY-потенциалом, которая удобна тем, что для нее интеграл (2) с фурье-образом (5) вычисляется точно.

Конструктивное сходство фурье-образов (5) и (6) позволяет выполнить необходимые преобразования и оценить вклад слагаемых сразу в двух моделях. Факторизация аргумента подынтегрального логарифма приводит к следующему его представлению:

где

Для сокращения записей переменная $x$, входящая в термодинамические величины, не снабжена индексами, поэтому следует иметь в виду, что в случае DY-потенциала она равна x = nβw_DY, а в случае DE-потенциала x = nβw_DE.

Подынтегральный логарифм преобразуется к виду

где y(t) = y₂(t)/y₁(t), y₁(t) = ln(f₁(t)), y₂(t) = |ln(1 – ‒ f₂(t)/f₁(t))|.

Поскольку Y(t) > 0, y₁(t) > 0, то из (7) следует, что 0 < y(t) < 1 при t > 0. Покажем, что функция y(t) монотонно убывает. Применению известных признаков препятствует трансцендентное уравнение y'(t) = 0, которое не допускает точного решения. Поэтому здесь монотонность y(t) устанавливается с помощью функции расстояния ρ(t) = = y₁(t) − y₂(t) между убывающими функциями y₁(t) и y₂(t). Функция ρ(t) достигает максимума в некоторой точке t₀ > 0. Это означает, что y(t) = 1 − − ρ(t)/y₁(t) может иметь минимум в точке t₀ или ее окрестности. Если последнее предположение верно, то при условии $\mathop {\lim }\limits_{t \to \infty } y\left( t \right) = 0$ следует ожидать существования максимума функции y(t) и соответствующего минимума функции ρ(t) при t > t₀. Однако ρ(t) имеет единственный экстремум – максимум в точке t₀. Следовательно, y(t) монотонно убывает при t → ∞, а ее значение

является наибольшим (графики функции y(t) при актуальных значениях параметров, используемых для расчетов свойств, представлены на рис. 2), и выполняется неравенство y(t) < y(0). О величине x = nβw, входящей в (8), до конкретных расчетов известно лишь, что x > 0 (малые значения x отвечают высоким температурам, большие x – низким). Поэтому представляет интерес оценка y(0), содержащая только параметры D и d, соответствующие некоторой точке (δ, ε) ∈ G₁. Такая оценка дается неравенством доказательство которого приведено в Приложении.

Таким образом, реализуется наиболее благоприятный из возможных сценариев поведения функции y(t), которая вносит основной вклад в Y(t) только в окрестности точки t = 0 (если бы y(t) имела несколько близких к единице максимумов при $t > 0$, то ситуация была бы менее обнадеживающей). Отбрасывая функцию y(t) в (7), приходим к приближенному равенству ln(1 + nβ${\tilde {v}}(t)$) ≈ ≈ ln(f₁(t)), позволяющему существенно упростить вычисления. Ответ на вопрос об эффективности предлагаемой аппроксимации дает сопоставление результатов приближенных расчетов с экспериментом, а для DY-модели, кроме того, с данными, полученными по “точным” формулам. Здесь и далее кавычки подчеркивают приближенный характер выражения (1) при точном вычислении интеграла (2) с фурье-образом (5), т.е. под “точными” следует понимать формулы DY-модели, получаемые после вычисления (2) с функцией Y(t) = ln(1 + nβ${\tilde {v}}(t)$). Под приближенными формулами понимаются формулы DY- и DE-моделей, найденные при интегрировании с функцией y₁(t) = ln(f₁(t)). В этом порядке ниже приводятся выражения для свободной энергии:

где ${{q}_{d}}\left( x \right) = \sqrt {1 + xd} ,$ ${{Q}_{{{\text{DE}}}}}\left( x \right) = \sqrt {1 + {{q}_{d}}\left( x \right)} {\kern 1pt} .$

Последовательность, в которой записаны формулы для F, сохранена далее и для других соотношений, в частности для уравнений состояния, получаемых с помощью известного соотношения P = −(∂F/∂V)_T:

Здесь

Расчеты свойств удобно проводить в приведенных переменных ω = n/n_c, τ = T/T_c, Π = P/P_c (индекс “с” указывает на принадлежность величин к критической точке).

Критическое значение x_c = n_cβ_cw, используемое для преобразования УС к безразмерной форме, является решением системы, определяющей критическое состояние

где

Система (12) сводится к нелинейному уравнению

которое получено путем исключения параметра L_c из уравнений системы (13) и в “точном” случае может быть решено только численно. Однако в приближенных случаях уравнение (13) принимает компактный вид:

– для DY-модели

– для DE-модели

и допускает аналитическое рассмотрение. Решение уравнения (14) можно найти приближенно в виде разложения по малому параметру ξ = 1 – ‒ (D/d) < 1. В первом порядке теории возмущений его решение имеет вид где

Приближенное решение (15) непригодно вблизи верхней границы ε = δ⁴ области G₁, где параметр D становится малым, а условие малости ξ нарушается.

Если значение x_c = x_c(δ, ε) установлено, то подставляя в (11) отношение L_c, найденное из первого уравнения системы (13), приходим к уравнениям состояния в безразмерной форме

где

Уравнения (16) не зависят от индивидуальных характеристик вещества, но содержат безразмерные параметры δ и ε.

Предсказательные возможности УС иллюстрируются сопоставлением результатов расчетов температурных зависимостей изобарной теплоемкости C_p(T), коэффициента Джоуля–Томсона α_P(T) и скорости звука u_P(T) при постоянном давлении с данными экспериментов:

Здесь

M – молярная масса вещества.

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

Химические потенциалы μ(ω, τ) определяются соотношениями

Функции q_D(x), q_1D(x), Q(x), Q₁(x), введенные после формул (10), (11), являются функциями ω и τ, так как x = x_cω/τ.

РАСЧЕТЫ И ЭКСПЕРИМЕНТ

Приближенный характер используемых соотношений не позволяет получить точное описание свойств реальных систем. Однако, учитывая, что соотношения (17)–(20) содержат настраиваемые параметры δ, ε, можно попытаться согласовать данные теории с экспериментом. Для этого необходимо определить точку (δ, ε) ∈ G₁, в которой отклонения расчетных значений от соответствующих экспериментальных данных окажутся минимальными.

В общем виде задача поиска оптимальной точки (δ, ε) ∈ G₁ формулируется как задача нахождения минимума функции двух переменных

где $X_{i}^{{\exp }}$ – экспериментальные и $X_{i}^{{{\text{theor}}}}$(δ, ε) – расчетные значения i-го свойства.

Из необходимых условий экстремума функции (∂Φ/∂δ = 0, ∂Φ/∂ε = 0) можно получить несколько систем нелинейных уравнений, вид которых зависит от входящих в Φ(δ, ε) свойств. При этом неизбежно возникают известные вопросы о существовании решений, о начальном приближении и др., требующие детального изучения в каждом случае. Однако здесь для нахождения минимума функции (21), по-видимому, целесообразно ограничиться непосредственным вычислением значений функции на дискретном множестве точек в секторе ${{G}_{1}}$, что позволит приближенно определить координаты оптимальных точек. Следует отметить, что глубина и координаты минимума поверхности Φ(δ, ε) зависят от набора свойств, привлекаемых для построения функции (21). В частности, если сумму Φ(δ, ε) образуют величины $u_{P}^{{{\text{min}}}}(\delta ,\varepsilon ),$ $\alpha _{P}^{{{\text{max}}}}(\delta ,\varepsilon ),$ найденные по приближенным УС с отвечающими им экспериментальными значениями $u_{P}^{{{\text{min}}}}$ = 257.0 м/с, $\alpha _{P}^{{{\text{max}}}}$ = 0.426 К/Па для аргона при P = 10 МПа [19, 20], то для DY-потенциала минимум Φ(δ, ε) достигается при δ_min= 0.7, ε_min= 0.288 (применяются в расчетах и по “точным” формулам для DY-потенциала), а для DE-потенциала — при δ_min= 0.2, ε_min= 8.34 × 10^–4 .

Значения плотности ω(τ), используемые для построения зависимостей (17)–(19) в температурном интервале от 100 до 300 К при P = 10 МПа определяются в каждом случае с помощью соответствующих УС с учетом координат оптимальных точек.

На рис. 3 показаны результаты расчета удельной изобарной теплоемкости C_p(T) при P = 10 МПа. Видно, что кривые 1 и 2, построенные по “точному” и приближенному УС для DY-потенциала, при малых температурах демонстрируют поведение, не соответствующее данным измерений – кривая 4, что связано с чрезвычайно резким, хотя и качественно верным возрастанием плотности ω(τ) при понижении температуры. Зависимость ω(τ) изображена на фрагменте рис. 3. Близкое расположение кривых 1 и 2 означает, что приближенное DY-УС хорошо воспроизводит особенности, характерные для “точного” DY-УС. Следовательно, выбор главной части ln(f₁(t)) подынтегральной функции при вычислении I(n, β) сделан правильно, и можно предположить, что воображаемая кривая C_p(T), которую можно было бы получить при наличии “точного” УС для DE-потенциала, будет близка к кривой 3, построенной по приближенному DE-УС и согласующейся с экспериментом существенно лучше, чем кривые 1, 2.

На рис. 4 изображены экспериментальная [19] и теоретические зависимости скорости звука от температуры при P = 10 МПа для аргона. При высоких температурах все три рассмотренных УС приводят к описанию данных измерений с погрешностями, почти не зависящими от температуры при Т = 300 К – δu_P(DY) ≈ 5%, δu_P(DE) ≈ 15%. С понижением температуры для всех обсуждаемых УС расхождения с экспериментом увеличиваются, а УС-DЕ показывает наихудший результат. Погрешности в окрестности точки минимума функции u_p(T) составляют δu_P(DY) ≈ 15%, δu_P(DE) ≈ 28%. В области низких температур все три УС дают увеличивающиеся с понижением температуры (и плотности) отклонения от экспериментальной кривой 4.

Сравнение результатов расчетов коэффициента Джоуля–Томсона (18), представленных на рис. 5, с данными эксперимента [20] показывает, что α_P(T) является единственной среди рассмотренных величин, которая удовлетворительно описывается всеми УС. Погрешности в окрестностях точек максимумов кривых α_P(T) составляют δα_P(DY) ≈ 18% и δα_P(DE) ≈ 11%.

На рис. 6 представлены результаты расчетов линий сосуществующих фаз и экспериментальная кривая [19]. Расчет газовой ветви кривой по DY-УС хорошо, по DE-УС превосходно описывает реальную картину, но для жидкостной ветви УС приводят к значительным погрешностям, возрастающим с понижением температуры. Важно отметить, что для жидкой фазы “точное” и приближенное DY-УС воспроизводят весьма близкие друг к другу результаты, которые уступают по точности соответствующим данным, полученным по DE-УС, но свидетельствуют о качестве приближенного DY-УС.

Серия расчетов, касающихся модели с DY-потенциалом, демонстрирует вполне удовлетворительное описание приближенным DY-УС результатов, полученных с помощью “точного” DY-УС, на что указывает взаимное расположение кривых 1 и 2 на рис. 3–6. Подобное сравнение для DE-модели оказывается невозможным, поскольку точное DE-УС не найдено. Однако возможна априорная оценка точности расчетов, выполненных по приближенному DE-УС. Для такой оценки удобно использовать функцию y(t), введенную после формулы (7), характеризующую работоспособность аппроксимации подынтегрального логарифма Y(t) = ln(1 + nβ${\tilde {v}}(t)$) функцией y₁(t) = = ln(f₁(t)), необходимой при вычислении интеграла I(n, β).

На рис. 2 приведены графики функции y(t), содержащей параметр $x = n\beta w,$ при значениях x = 10 и x = 0.5, которые соответствуют низким и высоким температурам, взятым в пределах исследованного интервала температур. Ход кривых 1, 2, построенных с учетом координат (δ, ε) оптимальных точек для каждой модели, полностью отвечает анализу поведения функции y(t), проведенному выше для рассматриваемых моделей. Из рис. 2 и формулы (7) следует, что вклад y(t) в Y(t) для DE-модели, по крайней мере, не больше, чем в случае DY-модели. Это означает, что приближенное равенство Y(t) ≈ ln(f₁(t)) для DE-модели не должно быть менее точным, чем для DY-модели. Следовательно, приближенное DE-УС предсказывает результаты, не уступающие по точности результатам расчетов по приближенному DY-УС относительно соответствующих “точных” данных.

ЗАКЛЮЧЕНИЕ

Основные результаты работы сводятся к следующему:

1. Получение “точного” УС системы с DE-потенциалом на основе интегрального представления свободной энергии F сопряжено с трудностями, связанными с вычислением соответствующего интеграла $I\left( {n,\beta } \right)$, и приводит к чрезвычайно громоздким соотношениям, неудобным для дальнейшего применения. Для упрощения задачи предлагается аппроксимация подынтегрального логарифма путем выделения главной части, удобной для интегрирования и сохраняющей свойства модельной системы.

2. Возможность использования аппроксимации, основанной на замене ln(1 + nβ${\tilde {v}}(t)$) на ln(f₁(t)) при вычислении I(n, β), следует из неравенства y(t) < 1, доказанного для DY- и DE-моделей. С целью иллюстрации эффективности такой замены привлекается DY-модель, для которой вычисления проведены по “точным” и приближенным формулам, передающим с хорошей точностью все особенности поведения модельной системы, возникающие при “точном” описании.

3. Поскольку для DE-модели “точное” УС неизвестно, то прямой расчет погрешностей в этом случае невозможен. Поэтому выполнена оценка точности приближенного равенства ln(1 + nβ${\tilde {v}}(t)$) ≈ ≈ ln(f₁(t)) для обеих моделей. С учетом представлений относительно поведения функции y(t), характеризующей точность приближенного равенства, и аналогии, проводимой между моделями, показано, что соотношения DE-модели сравнимы по точности с подобными соотношениями DY-модели. Это означает, что приближенное DE-УС воспроизводит результаты, близкие к тем, которые могли быть получены по “точному” DE-УС.

4. Сопоставление выполненных здесь расчетов свойств моделей с данными экспериментов показывает, что DE-УС демонстрирует физически разумные результаты, в отличие от DY-УС, которое приводит, в частности, к качественно неверной зависимости ${{C}_{P}}\left( T \right)$ при низких температурах. Поэтому DE-УС, полученное в рамках формализма, основанного на выражении для свободной энергии $F$, представляется более пригодным для описания равновесных термодинамических свойств простой жидкости.

5. Недоступными для исследований на основе рассмотренных УС являются области низких температур и высоких давлений, где существенным оказывается влияние многочастичных сил. Общий недостаток моделей состоит в неправильном описании поведения изохорной теплоемкости C_V(T).

6. Предложенный способ аппроксимации можно рекомендовать для нахождения УС систем с потенциалами взаимодействия, аналитическое исследование которых в рамках подхода Д.Н. Зубарева пока не представляется возможным.

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