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

ВВЕДЕНИЕ

Изучение фазовых диаграмм различных веществ и определение границ фаз являются одной из основных задач статистической механики и термодинамики [1, 2]. Соответствующие исследования продолжаются уже на протяжении более 100 лет. К настоящему времени для многих веществ линии фазовых границ уже достаточно точно измерены. Кроме этого, на основе строгих теорий разработан целый ряд методов численного моделирования. В рамках классической статистики для их применения нужно знать лишь межчастичный потенциал или потенциалы взаимодействия [3, 4]. Эти методы дополняются квантовыми методами моделирования, которые в наиболее продвинутых вариантах рассматривают вещество как кулоновскую смесь электронов и ядер, т.е. как оно и существует в природе [5]. Одной из первых работ, в которой подобный подход был реализован, является [6], а подробное описание самого метода можно найти в [7].

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

В экспериментах достаточно точно измеряются лишь температуры до ~5 кК [9], причем для однородной жидкой фазы. А в случае, если нужно определить линию фазового перехода жидкость–газ, то проблемы возникают и при более низких температурах, так как нужно изучать разные по плотности, но равновесные фазы одновременно. Здесь и далее ветвь бинодали с меньшей плотностью для краткости называется “газовой ветвью бинодали”, а ветвь с большей плотностью – “жидкой ветвью бинодали”. В результате из всех металлов надежные измерения по критическим точкам и бинодалям есть в настоящee время лишь для двух металлов – ртути и цезия, которые обладают сравнительно низкими критическими температурами. Для ртути это T_c = 1753 К [10], а для цезия T_c =1924–1938 К [11, 12].

Здесь и далее Т – температура, P – давление, ρ – массовая плотность, n – концентрация. Индекс “с” относится к величинам в критической точке. Для других же щелочных металлов T_c расположена выше 2000 К. Поэтому для них оценки координат критической точки имеют больший разброс даже при наличии измерений положения бинодали [13, 14]. Для всех же остальных металлов прямых измерений положения бинодалей нет. Это относится и к ряду кристаллических (при нормальных условиях) диэлектриков и полупроводников. Тем не менее и для таких веществ возможны оценки при сравнительно низких температурах. В области газа имеются справочные данные по давлению насыщенных паров P(T) при давлениях до 1 атм [15, 16]. При этих условиях металл (или полупроводник) с высокой точностью является идеальным атомарным газом, что позволяет определить его плотность по уравнению состояния идеального газа. Это дает оценку ρ_Gas(T) на газовой ветви бинодали. В области жидкости для металлов есть измерения плотности на изобарах [9]. Так как жидкость обладает очень малой сжимаемостью, то эти плотности практически располагаются на жидкой ветви бинодали. Причем чем выше давление, тем менее заметно различие между бинодалью и изобарой (см., например, [10]). Это дает достаточно надежную оценку плотности на жидкой ветви бинодали ρ_Liq(T). Но эти оценки плотностей на обеих ветвях бинодалей ограничиваются, как правило, сравнительно низкими температурами, обычно порядка температуры кипения (при 1 атм), которые в 2–3 раза ниже, чем предполагаемая критическая температура [16]. Например, для того же цезия температура кипения – 670 К, а для ртути – 630 К.

Здесь на помощь, казалось бы, должны были прийти теория и расчеты. И действительно, для жидких металлов разработано достаточно много моделей взаимодействия. Для моделирования в рамках классической статистики наиболее продвинутой моделью взаимодействия является метод погруженного атома (Embedded Atom Me-thod) [17–19]. Далее для него используется общепринятое сокращение EAM. В рамках EAM исследуемое вещество представляется как однокомпонентная система, состоящая из точечных беcструктурных частиц. Для нее строится многочастичный потенциал взаимодействия, параметры которого выбираются так, чтобы описывать конкретное вещество. Впервые такая процедура была реализована в работах [20, 21]. Далее методами Монте-Карло или молекулярной динамики можно рассчитать различные термодинамические и переносные свойства такой системы. Так, в частности, были рассчитаны и бинодали ряда металлов [19]. Однако бинодали именно модельной системы, но не реального вещества. Это связано с тем, что, хотя область температур, описываемых EAM, может быть существенно выше температуры кипения изучаемого вещества, этот потенциал неприменим для газовой фазы [17, 18]. В газовой же фазе различных металлов при не слишком больших плотностях хорошо работают парные потенциалы для димеров [22]. Но, очевидно, они неприменимы при увеличении плотности и в жидкой фазе. Таким образом, для металлов не существует единого потенциала взаимодействия, применимого во всех фазах, что делает невозможным прямой расчет в рамках классической статистики. Данная проблема могла быть решена в рамках квантовых вычислений, но они тоже хорошо описывают, как правило, достаточно плотные системы [5–7]. Помимо этого, для них существуют и фундаментальные проблемы, связанные, например, с неопределенностью обменно-корреляционного взаимодействия, которые только усугубляются при понижении плотности. Также они требуют значительных вычислительных ресурсов, поэтому число частиц, используемое сейчас в квантовoм моделировании, может быть также недостаточным для моделирования двух фаз. В результате в рамках этих методов в настоящий момент можно рассчитать критическую изотерму и оценить координаты критической точки [23, 24], но прямой расчет бинодали все еще невозможен.

Таким образом, для металлов напрямую получить бинодаль из экспериментов или расчетов пока невозможно. Поэтому приходится использовать менее точные методы, построенные на соотношениях подобия и экстраполяции известных из эксперимента низкотемпературных данных. Этот подход развивается уже несколько десятилетий, и подобным образом были оценены критические точки для многих веществ, включая металлы [25–27]. Одними из наиболее общих соотношений подобия являются соотношения, связанные с линией единичного фактора сжимаемости на плоскости плотность–температура [28–32]. На их основе ранее было построено общее уравнение для бинодали [33, 34], которое было успешно проверено для веществ с известной бинодалью [34]. Также этот подход уже применялся для ряда металлов, включая алюминий, медь, уран и галлий [35, 36]. Цель данной работы – используя данный метод, построить бинодаль висмута и оценить положение его критической точки.

ЛИНИЯ ЕДИНИЧНОГО ФАКТОРА СЖИМАЕМОСТИ И БИНОДАЛЬ

Линия единичного фактора сжимаемости (Zeno line), как следует из ее названия, это линия, на которой сжимаемость системы Z = 1. Так как Z = = P/(nk_BT), где k_B – постоянная Больцмана, то на этой линии давление любой системы будет совпадать с давлением идеального газа при одних и тех же n и T. Это означает, что вдоль данной линии силы притяжения и силы отталкивания взаимно компенсируются. Очевидно, что для систем с чисто отталкивательным взаимодействием ее не существует. Как оказалось, она может обладать одним очень интересным и универсальным свойством, а именно, эта линия для очень многих модельных и реальных систем является прямой во всей области существования жидкости от газовых плотностей до линии плавления [32]. Впервые это прямолинейное поведение было обнаружено Бачинским, когда он применил уравнение Ван-дер-Ваальса (УВ) для описания своих экспериментов [37]. Поэтому сначала считалось, что такое поведение является свойством только этого уравнения. И о ней забыли на несколько десятилетий. Но потом анализ экспериментальных данных для целого ряда реальных веществ, которые не описываются УВ, подтвердил, что контур Z = 1 является прямой линией и для них. Например, в базе данных [8] только у пяти веществ эта линия – непрямая. Для определения прямой линии достаточно двух параметров, а именно координат точек пересечения с осями температуры и плотности. И ее уравнение записывается как

Для УВ Т_B = a/b, n_B = 1/b, где a, b – константы уравнения Ван-дер-Ваальса. Для других же систем, очевидно, требуются дополнительные методы. Самым простым и надежным теоретическим методом является вириальное разложение [1, 2]. Оно гарантировано работает практически для любой системы при малых плотностях. С его помощью можно также получить прямую линию для контура Z = 1 [29, 32]. При этом параметры Т_B, n_B выражаются через второй и третий вириальные коэффициенты, причем Т_B оказывается не чем иным, как температурой Бойля. Поэтому величины Т_B, n_B получили название бойлевских параметров. И хотя для высоких плотностей вириальное разложение может быть уже неприменимо, линия Z = 1, получаемая в численных расчетах для модельных систем, в точности совпадает с линией, получаемой из вириального разложения [32]. Отметим также, что линия Z = 1 является лишь одной из идеальных линий, т.е. таких, на которых некоторое свойство вещества совпадаeт с аналогичным свойством идеального газа. Поэтому существуют идеальные линии, например, для энтальпии и внутренней энергии, причем первая тоже прямая [32, 38]. Заметим также, что идеальные линии можно рассмотреть и в кристаллической фазе, но там они уже не являются прямыми [39].

Еще одним свойством линии Z = 1 является ее положение относительно бинодали. Из общих термодинамических соотношений было показано, что линия Z = 1 является касательной к жидкой ветви бинодали при Т → 0 [28, 32]. Причем это свойство не зависит от того, является линия Z = 1 прямой или нет. Общий вид фазовой диаграммы жидкости c линией Z = 1 в координатах, приведенных к критическим, показан на рис. 1 на примере системы Ван-дер-Ваальса. И качественно аналогичную фазовую диаграмму можно наблюдать у многих других систем [29, 32].

Таким образом, линия Z = 1 является асимптотой жидкой ветви бинодали при низких температурах. Это дает линейную зависимость плотности на этой ветви бинодали от температуры при Т → 0. Для газовой же ветви в том же пределе, согласно уравнению Клапейрона–Клаузиуса, получается экспоненциальная зависимость плотности от температуры [1, 2]. В другом пределе, при Т → Т_с, поведение обеих ветвей следует из теории скeйлинга, т.е. это степенная зависимость с нецелым показателем степени [1]. Кроме этого, можно использовать симметризацию бинодали относительно ее диаметра на основе идеи изоморфизма, предложенной в [30]. Все эти свойства позволили построить общее уравнение для обеих ветвей бинодали [34], которое имеет вид

В (2) индексы L и G обозначают плотность на жидкой и газовой ветвях бинодали соответственно. Индекс 2D обозначает удвоенную плотность на диаметре бинодали, т.е. n_2D = n_L + n_G. Асимптотическое поведение при Т → 0 определяет выражение для коэффициентов A и В. А показатель степени β = 0.326 – критический индекс, который определяется теорией скeйлинга [1]. Величина же q = Q/(k_BT_с), где Q – это некоторая эффективная теплота испарения. Величина q изначально рассматривалась как подгоночный параметр. В [34] этот параметр определялся из условия наименьшего отклонения плотностей, рассчитанных по (2), от экспериментальных плотностей для веществ с известными бинодалями из базы данных [8]. Ошибка такого сравнения была менее 4%. Помимо этого, обнаружилось, что Q имеет значения, очень близкие к теплоте испарения при P = 1 атм. Таким образом, проверка уравнения (2) показала, что оно может успешно описать зависимость плотностей на обеих ветвях бинодали для широкого класса веществ. Поэтому возникает идея применить его для веществ, для которых бинодаль неизвестна, например для металлов. Но в этом случае в (2) входят пять неизвестных параметров, которые должны быть определены дополнительно. Это два бойлевских параметра, плотность и температура в критической точке и q. Об определении первых двух будет сказано ниже. Что же касается оставшихся, то их число может быть уменьшено, так как ранее [32–34, 40] было найдено соотношение, связывающее бойлевские и критические параметры:

где S₁ ≈ 0.67. Соотношение (3) выполняется практически для всех веществ из [8] (см. табл. 1 в [40]), поэтому считается возможным применить его и в данном исследовании. Если T_B и n_B известны, то в (2) остается только два параметра: q и T_с или q и n_с. Для их определения в [34] была разработана соответствующая процедура минимизации. Для нее, однако, нужны некоторые данные на части бинодали. О выборе этих данных и о процедуре минимизации будет сказано ниже. Сначала следует отметить еще два важных пункта.

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

Здесь Z_с – фактор сжимаемости в критической точке, соответственно, давление P_с = (n_с)²k_BT_с/n_B. Уравнение (4) известно как соотношение Тиммерманса [28, 41]. Оно имеет несколько меньшую область применимости, чем уравнение (3). Тем не менее во многих случаях оно дает достаточно точные значения критического давления. Поэтому будем использовать его и для оценки критического давления висмута.

Вторым важным пунктом является то, что любой металл в жидком состоянии является двухкомпонентной системой, состоящeй из электронов и положительных ионов. Это может приводить к неоднозначному определению самой величины фактора сжимаемости Z [32]. Но если рассматривать металл как однокомпонентную систему с некоторым эффективным потенциалом, то можно пользоваться обычным определением Z = = P/(nT) = m_aP/(ρT), m_a – масса частицы (или атома). Это определение и используется далее.

ПРИМЕНЕНИЕ СООТНОШЕНИЙ ПОДОБИЯ ДЛЯ ВИСМУТА И ЕГО БИНОДАЛЬ

Висмут обладает очень низкой температурой плавления при нормальных условиях – Т_m = 544.7 К, которая приблизительно соответствует его тройной точке. Его плотность при нормальных условиях ρ₀ = 9.79 г/см³. При этом он обладает “обратным” плавлением, т.е. плотность жидкости в точке плавления больше, чем плотность твердой фазы ρ₀, и составляет 10.05 г/см³. Температура кипения висмута – 1837 К при P = 1 атм [16]. Температура сравнительно низкая, поэтому в области жидкости до температур порядка температуры кипения существует достаточно много экспериментальных данных, например [42]. Эти данные позволяют построить достаточно точное аналитическое уравнение состояния для давления. Поэтому линия Z = 1 определена в настоящей работе непосредственно по данным [42], но лишь до Т = 973 К и P = 4.3 ГПа. И она оказалась прямой в этой области. При более высоких температурах эту линию можно рассчитать, используя EAM, как сделано для галлия в [36]. Соответствующий потенциал представлен в [43]. При Т < 973 К линия Z = 1, рассчитанная по EAM, совпадает с данными, полученными из экспериментов [42]. При высоких температурах расчетная линия также является прямой, дополняющей данные при Т < 973 К. Однако, как было сказано во Введении, расчеты по ЕАМ тоже ограничены по температуре сверху, поэтому здесь решено ограничиться температурами ниже 3000 К.

Для расчетов использовалось моделирование в рамках метода Монте-Карло в каноническом ансамбле [3]. Этот метод позволяет рассчитать давление как функцию плотности и температуры. Поэтому, варьируя плотность на различных изотермах при T ≤ 3000 К, были найдены точки, где Z = 1. В настоящем моделировании использовалось N = 1000 частиц, случайно размещенных в кубической ячейке c длиной ребра L = (N/n)^1/3. Для достижения равновесия использовались 300 конфигураций (одна конфигурация соответствует смещению или несмещению всех N частиц). Еще 700 конфигураций использовались для осреднения. Радиус обрезания полагался равным L/2. За его пределами использовано стандартное условие равенства парной корреляционной функции единице [3], необходимое для парной части потенциала.

В рассматриваемой области температур и плотностей к полученным результатам применен метод наименьших квадратов и определено, что построенная таким образом линия Z = 1 является прямой с точностью ~1%. Для определения бойлевских параметров полученная зависимость экстраполирована к нулевым значениям плотности и температуры и найдено, что для висмута T_B = = 11 510 К, ρ_B = m_an_B = 10.636 г/см³ (m_a – масса атома висмута).

Таким образом, можно применить упомянутую выше процедуру минимизации, но сначала нужно определиться с экспериментальными данными, относительно которых она будет проводиться. Таковыми, очевидно, являются данные измерений при низких температурах (см. Введение). Давление насыщенного пара висмута представлено в [16] до Т = 1695 К, что близко к точке кипения при нормальных условиях. Используя уравнение идеального газа, из этих данных можно получить плотность на газовой ветви бинодали. Зависимость плотности жидкости от температуры при P ~ 1 атм измерена в [44] до Т = 1600 К. Здесь она принимается за плотность на жидкой ветви бинодали. Далее значение плотностей на обеих ветвях бинодали вычисляется по (2) с учетом (3) до температуры 1600 К при некоторых q и T_с. Варьируя эти два параметра, можно найти их оптимальные значения, при которых разница между плотностями по (2) и экспериментальными данными является наименьшей. Соответствующие формулы для q и T_с приведены в [34]. Теперь можно построить всю бинодаль от тройной до критической точки. Критическое давление было определено по соотношению Тиммерманса (4). Величина q = 6.2, что дает Q = 164 кДж/моль, а справочная теплота парообразования при P = 1 атм – 177 кДж/моль [16]. Остальные критические параметры представлены в таблице наряду с оценками других авторов.

Как можно видеть из таблицы, T_c данной работы ближе к данным [26], а Z_c – к данным [27]. При этом настоящие критические плотность и давление – наибольшие из всех оценок. Для щелочных металлов Z_c ~ 0.2, а для разнообразных газов Z_c ~ 0.3, еще большее значение Z_c ~ 0.36 характерно для ртути [32]. Бинодаль вместе с линией Z = 1 представлена на рис. 2. Туда же добавлена кривая плавления согласно данным измерений [45].

В [45] представлены данные по давлению, температуре и плотности твердой фазы при плавлении, но нет плотности жидкой фазы, кроме как при Т = Т_m. Для определения последней взяты те же температура и давление, что и в [45], и рассчитаны соответствующая плотность в жидкой фазе с использованием потенциала [43] (рис. 2). При Т = = Т_m полученная в расчете плотность жидкой фазы 10.07 г/см³ находится в хорошем согласии с экспериментальным значением, приведенным в [45], – 10.05 г/см³. Также отлично виден “обратный ход” плавления, при котором плотность твердой фазы меньше, чем плотности жидкости при одной и той же температуре. Он продолжается, пока температура не понизится до 461 К, а потом начинается обычный ход плавления.

ЗАКЛЮЧЕНИЕ

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

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