Журнал неорганической химии, 2023, T. 68, № 2, стр. 191-202

ВВЕДЕНИЕ

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

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

Температурные зависимости объема при атмосферном давлении могут быть описаны через коэффициент термического расширения в виде полиномиальных функций [1]. Термическое расширение также может быть представлено полуэмпирическими квазигармоническими моделями в диапазоне от 0 K до температуры плавления [2]. При параметризации УС часто используют информацию об объеме при различных температурах и давлениях. В литературе представлены зависимости V(T) для большого количества соединений. Высокая точность определения объема и температуры при атмосферном давлении обусловливает хорошее качество предложенных моделей и расширяет набор веществ, для которых предложены УС. Данные по изотермической сжимаемости описаны в литературе с помощью различных изотермических УС; одним из наиболее популярных является уравнение Берча–Мурнагана третьего порядка, которое до сих пор используется для описания P–V-соотношений веществ различных типов [3–6]. Наиболее востребованные изотермические УС были проанализированы в работах [7, 8]. Показано, что экспериментальные данные по изотермической сжимаемости могут быть описаны в пределах погрешности измерения до очень высоких давлений с помощью трех или четырех параметров, один из которых может быть зафиксирован. Для изотермических условий в работах [7, 9] рекомендовали применять модифицированное уравнение Тайта [10], которое использовано и в настоящей работе.

В то время как для описания изотермических данных предложено много различных УС, для моделирования свойств в широком интервале температур и давлений используется существенно меньшее число аналитических зависимостей [11]. Многие авторы для получения PVT-уравнений добавляют к изотермическому УС термическое давление [12–15], которое можно выразить различными способами: от простой линейной зависимости до сложного интеграла, например, как в уравнении Ми–Грюнайзена–Дебая [16].

Простейшим способом описания теплоемкости является полиномиальная зависимость от температуры. В литературе представлены различные аналитические зависимости, применяемые непосредственно (например, [17]) или в виде линейных комбинаций [18]. Помимо полиномиальных функций для описания температурной зависимости теплоемкости могут быть использованы функции Эйнштейна и Дебая [19]. Достаточно распространенным способом моделирования теплоемкости является подход CALPHAD третьего поколения, в котором теплоемкость при фиксированном давлении представляется как комбинация эйнштейновской модели с полиномиальной коррекцией на ангармоничность и наличие электронных эффектов [20]. В работе [21] предложен метод описания температурной зависимости теплоемкости с помощью комбинации функций Планка–Эйнштейна и показано, что этот подход позволяет описывать весь набор экспериментальных данных в пределах погрешности их определения, при этом количество параметров модели обычно не превышает 9–10. Такой метод был применен для аппроксимации термохимических измерений различных веществ [20, 22, 23]; при этом особо подчеркивалось, что параметры модели не имеют строгого физического смысла, однако возможность адекватного описания результатов измерений и корректное предельное поведение свойств делают его перспективным для использования в термодинамических базах данных.

Практический интерес представляет описание разнообразных наборов термодинамических данных единым УС. В литературе представлено несколько работ, в которых реализован такой подход. Авторы [13] определили параметры ранее упомянутого уравнения Ми–Грюнайзена–Дебая с помощью как PVT-данных, так и сведений по адиабатическому модулю упругости (K_S). Однако другие термодинамические свойства (теплоемкость, приращения энтальпии) не были включены в оптимизацию. Даже с учетом возможности расчета этих свойств качество их описания такими уравнениями оставляет желать лучшего. Аналитические зависимости сложны и при этом не позволяют воспроизвести термохимические данные с экспериментальной точностью (доли процента).

Если при построении УС исходить из аналитических зависимостей энергии Гиббса [24, 25] или Гельмгольца [26, 27], то с помощью подобных уравнений можно описать как PVT-соотношения, так и любые термодинамические свойства веществ. Например, авторы [28] провели параметризацию предложенного ими полуэмпирического уравнения для энергии Гельмгольца и получили для ряда веществ хорошее описание значений теплоемкости, коэффициента термического расширения, адиабатического модуля упругости при атмосферном давлении, а также предсказали с приемлемой точностью имеющиеся PVT-данные при высоких давлении и температуре.

Разрабатывать УС на основе комбинации эйнштейновских функций предлагалось ранее в работах [24, 27]. Целью настоящей работы является анализ возможности расширения метода Воронина–Куценка [21] для построения УС кристаллической фазы. Помимо изобарной теплоемкости с помощью функций Планка–Эйнштейна можно описывать также и изохорную теплоемкость [21]. Поэтому в настоящей работе приведены два варианта УС на основе комбинации функций Планка–Эйнштейна, в основу которых положены зависимости энергии Гиббса или энергии Гельмгольца от естественных переменных, а также проведено сравнение интерполяционных и экстраполяционных характеристик предложенных УС.

В качестве тестовой системы был выбран оксид магния (периклаз), для которого в литературе представлено большое количество различных экспериментальных данных в широком диапазоне давлений и температур. В процессе оптимизации уравнений состояния использовали следующие термодинамические свойства: изобарную теплоемкость и приращение энтальпии, термическое расширение, PVT-соотношения и адиабатический модуль упругости при различных температурах и давлениях.

РАСЧЕТНАЯ ЧАСТЬ

Как было показано ранее [21], применение комбинаций функций Планка–Эйнштейна демонстрирует удовлетворительные результаты в описании функций C_P(T), $S_{P}^{^\circ }(T)$ и H°(T) – H°(T_ref) кристаллических соединений. Расширить этот подход для описания различных свойств кристаллических веществ при изменении давления (или объема) можно с помощью основных соотношений химической термодинамики. В данной работе рассмотрены два способа: первый основан на параметризации функции C_P(P, T), которую можно преобразовать в зависимость энергии Гиббса от давления и температуры; второй – на параметризации функции ${{C}_{V}}$(V,T), которую можно преобразовать в зависимость энергии Гельмгольца от объема и температуры. Далее эти модели будут обозначаться как C_P-подход и ${{C}_{V}}$-подход, соответственно. Взаимосвязь между термодинамическими функциями и последовательность математических преобразований схематически представлены на рис. 1.

C_P-подход. Исходное выражение для изобарной теплоемкости C_P в рамках модели Воронина–Куценка [21] имеет вид:

где $x_{i}^{P}$ = ${{\left( {\theta _{i}^{P}} \right)} \mathord{\left/ {\vphantom {{\left( {\theta _{i}^{P}} \right)} T}} \right. \kern-0em} T};$ $\alpha _{i}^{P}$ и $\theta _{i}^{P}$ – варьируемые параметры, надстрочный символ P показывает отнесение параметра к C_P-варианту модели. Соответственно, энтропия может быть рассчитана по уравнению (2):

С помощью одного из соотношений Максвелла можно связать энтропию и объемные свойства:

Далее, используя выражение для S(T) в виде линейной комбинации функций Планка–Эйнштейна, можно описать термическое расширение, интегрируя соотношение Максвелла:

где V(P, T) – V(P, 0) – это изменение объема при постоянном давлении, V(P,0) представляет собой соответствующий уровень отсчета. Зависимость энтропии от давления в выражении для S_P(T) можно учесть, например, через зависимость от давления параметров $\alpha _{i}^{P}$ или $\theta _{i}^{P}.$ Таким образом, уравнение (2) преобразуется в обобщенную форму уравнения состояния Планка–Эйнштейна (уравнение (5)), которое далее будет обозначаться как ΔV_PE(P,T):

Анализ предложенного УС показал, что параметр $\alpha _{i}^{P}$ целесообразно принять независимым от давления, т.е. ${{\left( {\frac{{\partial \alpha _{i}^{P}}}{{\partial P}}} \right)}_{T}} = 0.$

Для определения других термодинамических свойств необходимо воспользоваться производными V(P, T) по давлению и температуре:

Изохорная теплоемкость ${{C}_{V}}$, изотермический (K_T) и адиабатический (K_S) модули упругости могут быть рассчитаны с помощью соотношений (8)–(10):

${{C}_{V}}$-подход. В качестве исходной для последующих преобразований может быть выбрана не только изобарная, но и изохорная теплоемкость:

где $x_{i}^{V}$ = ${{\left( {\theta _{i}^{V}} \right)} \mathord{\left/ {\vphantom {{\left( {\theta _{i}^{V}} \right)} T}} \right. \kern-0em} T};$ $\alpha _{i}^{V}$ и $\theta _{i}^{V}$ – варьируемые параметры, надстрочный символ V показывает отнесение параметра к ${{C}_{V}}$-варианту модели. Следовательно, зависимость энтропии от температуры при фиксированном объеме описывается выражением:

Если воспользоваться соотношением Максвелла:

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

Необходимые для последующих преобразований частные производные давления по объему и температуре записаны ниже (уравнения (15), (16)):

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

и приращение энтальпии как интеграл от функции C_P(Т). При расчете адиабатического модуля упругости K_S использованы те же уравнения (уравнения (9), (10)), что и в C_P-подходе. В настоящей работе эти расчеты осуществлялись с помощью численных методов.

Гибкость модели. В записанных выше интегральных уравнениях присутствуют функции V(P, 0) и P(V, 0), представляющие собой уровни отсчета. Поскольку большинство экспериментальных данных определено при комнатной температуре или выше, за референсную температуру удобнее выбрать 300 K; в этом случае требуется провести элементарные алгебраические преобразования:

где ΔV_PE(P, T) и ΔP_PE(V, T) рассчитываются по уравнениям (5) и (14).

Слагаемые V(P, 300) и P(V,300) в уравнениях (18) и (19) могут быть описаны любым изотермическим уравнением состояния. В настоящей работе для этого было использовано уравнение Тайта, модифицированное авторами [10]. Основными достоинствами этого уравнения являются аналитическая инвертируемость и хорошие аппроксимационные характеристики [7]:

В записанных выше уравнениях V₀ – объем при давлении P = 1 бар и температуре $T = 300\,\,{\text{K;}}$ K₀, $K_{0}^{'},$ – изотермический модуль упругости и его первая и вторая производные соответственно. При равенстве нулю параметра уравнение Хуан и Чоу преобразуется в хорошо известное уравнение Мурнагана [29]:

Варьируемым параметром предложенной нами модели является параметр θ_i, зависящий от давления (С_P-вариант) или объема (${{C}_{V}}$-вариант). В первом случае:

где B_i, C_i – дополнительные варьируемые параметры, $\theta _{i}^{{0P}}$ соответствует θ_i в изначальном варианте модели Воронина–Куценка [21], P₀ = 1 бар.

Во втором случае вид зависимости $\theta _{i}^{V}(V)$ был аналогичен использованному в работе [30]:

где $\gamma _{i}^{0}$ и q_i – варьируемые параметры, V₀ – объем при давлении P = 1 бар и температуре T = 300 K.

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

Оптимизация параметров модели. Два варианта уравнения состояния, предложенные в настоящей работе, включают следующие варьируемые параметры: V₀, K₀, $K_{0}^{'},$ для изотермической части, несколько пар α_i, $\theta _{i}^{0},$ описывающих теплоемкость при постоянном давлении (P = 1 бар) или постоянном объеме (объем при давлении P = = 1 бар и выбранной температуре), и несколько дополнительных параметров для учета зависимости θ_i от давления (B_i, C_i) или объема ($\gamma _{i}^{0},$ q_i) в зависимости от использованного подхода. Для получения достоверных оценок параметров необходимо проводить оптимизацию с использованием как можно большего числа разных типов экспериментальных данных. В настоящей работе при оптимизации учитывали результаты измерений теплоемкости и приращения энтальпии при P = 1 бар, объемные свойства и адиабатический модуль упругости при различных давлениях и температурах.

При параметризации обоих вариантов УС использовали идентичную целевую функцию, которая представляла собой сумму квадратов отклонений экспериментальных значений от расчетных с учетом статистических весов для каждого типа данных. Для C_P-подхода параметры оценивали следующим образом: при аппроксимации теплоемкости и приращения энтальпии рассчитывали параметры α_i и $\theta _{i}^{0},$ далее по изотермическим данным при 300 K оценивали параметры V₀, K₀, $K_{0}^{'},$ после чего находили полный набор параметров с учетом всех доступных экспериментальных данных, используя ранее подобранные значения параметров в качестве начального приближения. В случае ${{C}_{V}}$-варианта модели предварительно можно было оценить только значения параметров изотермического УС, все остальные параметры рассчитывали при оптимизации всего набора экспериментальных данных. Параметр V₀ был изначально фиксирован для обоих вариантов модели и принят равным 74.71 ų как результат усреднения значений молярного объема из разных источников [12, 14, 31, 32].

Рассчитанные параметры охарактеризованы 95%-ными доверительными интервалами; все рекомендованные значения являются статистически значимыми. Дополнительные параметры, не представленные в финальном наборе, принимали равными B_i = 0 или C_i = 1 (C_P-вариант), $\gamma _{i}^{0}$ = 0 или q_i = 1 (${{C}_{V}}$-вариант). Для адекватного описания экспериментальных данных оказалось достаточно четырех–пяти пар параметров α_i и θ_i. Модель с тремя парами недостаточно хорошо описывала теплоемкость, а при использовании более шести пар параметров некоторые из них становились статистически незначимыми. Было также замечено, что дополнительные параметры B_i, C_i и $\gamma _{i}^{0},$ q_i требуются в первую очередь в термах с наибольшими значениями $\theta _{i}^{0}.$ Введение зависимости от давления или объема в термы с наименьшими значениями $\theta _{i}^{0}$ не влияло на качество описания в целом, так как эти параметры “отвечают” за низкотемпературные вклады, а экспериментальные данные при высоком давлении и низкой температуре отсутствуют.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Численные значения параметров для C_P- и ${{C}_{V}}$-вариантов УС приведены в табл. 1; в первом случае для адекватного описания различных типов данных потребовалось 19, а во втором – 13 параметров. Общая сводка информации об источниках экспериментальных значений и точности аппроксимации данных предложенными моделями представлена в табл. 2.

Визуально о качестве предложенных моделей можно судить по рис. 2–8. На рис. 2 представлены экспериментальные данные по теплоемкости при атмосферном давлении; символами обозначены измеренные значения, линиями – результаты расчета. Видно, что низкотемпературные данные [33], полученные методом адиабатической калориметрии, вполне согласуются с результатами ДСК-опытов [34, 35] и значениями приращения энтальпии, полученными с помощью калориметрии сброса [36–38] (рис. 2, 3). Выше 100 K предложенные УС описывают теплоемкость с точностью до 2%, а приращение энтальпии – до 1%. Если в случае энтальпии точность совпадает с погрешностью экспериментального определения, то для теплоемкости разница между экспериментальными и расчетными данными выше, чем заявлено в оригинальных работах [34, 35] (0.5% для интервала 10 < T < 20 K и 0.2% для T > > 20 K). Однако следует отметить, что разброс значений C_P, полученных разными авторами, превышает указанные погрешности. Учитывая, что в качестве погрешности обычно приводят разницу между измеренными и рекомендуемыми значениями теплоемкости стандартов, полученное расхождение экспериментальных и расчетных значений можно считать удовлетворительным.

Расхождение между измеренными и сглаженными значениями C_P(T) при низких температурах может быть уменьшено за счет увеличения количества параметров $\alpha _{i}^{0}$ и $\theta _{i}^{0}.$ В настоящей работе мы ограничились представленными наборами, так как ошибка 3–4% в описании теплоемкости при T < 50 K оказывает незначительное влияние на значения стандартных термодинамических функций при комнатной температуре и выше [39].

Оба варианта УС показывают близкое описание данных по термическому расширению периклаза до 3000 K; расхождение между расчетом и экспериментом не превышает ошибку эксперимента, соответствующую 0.5% [1]. На рис. 4 изображены результаты измерений объема MgO в диапазоне температур 300–3000 K при стандартном давлении P = 1 бар; видно хорошее согласие между данными разных авторов [1, 40] и модельным описанием. Результаты работ [12, 31, 32, 41–43] по изотермической сжимаемости (300 K) MgO также хорошо согласуются между собой (рис. 5) и с результатами расчетов с помощью предложенных вариантов УС. В случае изучения объемных свойств при высоких давлениях есть определенные проблемы, связанные непосредственно со способом определения значения P; его рассчитывают на основании свойств вещества сравнения (так называемого калибранта давления). В работе [43] показано, что использование разных уравнений состояния одного и того же вещества сравнения при фиксированной температуре приводит к разнице в 20 ГПа при абсолютном значении 140 ГПа. С учетом этого разница рассчитанных и измеренных PVT-данных на рис. 6 не представляется существенной, она соответствует возможной разнице в шкалах веществ сравнения [43]. Следует также подчеркнуть, что C_V-вариант УС лучше описывает PVT-данные при высоких давлениях.

С высокой точностью объемные свойства можно определить через адиабатический модуль упругости, так как измеряемые значения не зависят от веществ сравнения. На рис. 7а представлена температурная зависимость адиабатического модуля упругости K_S при атмосферном давлении. Оба подхода (C_P и C_V) хорошо описывают результаты измерений [44–47], укладываясь в погрешность эксперимента 1% согласно данным [47] (рис. 7б). Адиабатический модуль упругости при высоких давлениях можно получить комбинацией сверхзвуковой интерферометрии и рентгеновской дифракции; в работах [14, 48] такие измерения были выполнены до температуры 1650 K и давления 23 ГПа. На рис. 8 представлены относительные отклонения измеренных и рассчитанных значений K_S при давлениях до 24 ГПа и температуре от 300 до 1650 K. Оба варианта УС демонстрируют близкое описание K_S–V–T данных; относительные отклонения не превышают 2%, что немного превышает экспериментальную погрешность 1.2%, указанную в работе [14].

Таким образом, можно сделать вывод, что обе модели демонстрируют в целом схожее и адекватное описание эксперимента, в то время как C_P-подход немного точнее описывает данные при атмосферном давлении, а C_V-подход – данные при высоком давлении.

Предсказательная способность. Стандартным способом проверки предсказательной способности любой модели является воспроизведение данных, не использованных при ее параметризации. В настоящей работе для такой проверки при оптимизации были: а) исключены из общей выборки данные по адиабатическому модулю упругости при высоких давлениях, б) учтены результаты измерений коэффициентов термического расширения, адиабатического модуля упругости, теплоемкости и приращения энтальпии при стандартном давлении, а также объемные свойства при 300 K, в) учтены только изотермические PV-данные, теплоемкости, приращения энтальпии и значения коэффициента термического расширения при P = 1 бар. После получения статистически значимых параметров обеих моделей были рассчитаны средние относительные отклонения для разных типов литературных данных.

Вариант (а). В табл. 2 приведены характеристики модельного описания для полного (Опт-1) и ограниченного (Опт-2) наборов данных. Видно (столбцы Опт-1(C_P) и Опт-2(C_P)), что в случае C_P-варианта УС, построенного по усеченному набору данных, получается менее точная оценка значений K_S при высоком давлении, в то время как другие свойства воспроизводятся более точно. Частичное ухудшение описания теплоемкости и энтальпии связано с меньшим количеством использованных пар параметров α_i и $\theta _{i}^{0}.$ Параметры УС, рассчитанные при оптимизации ограниченного набора данных (Опт-2), перечислены в табл. S2 . При использовании C_V-варианта модели значения MRD практически одинаковы при полном и усеченном наборе данных (столбцы Опт-1(C_V) и Опт-2(C_V)).

Вариант (б). В этом случае (Опт-3) наблюдалось приемлемое описание PVT-свойств при высоких температурах и давлениях. Данные по изотермической сжимаемости были использованы для получения параметров изотермического УС, которые затем фиксировали при последующей оптимизации. Лучшие варианты описания (с 4 или 5-ю парами α_i и $\theta _{i}^{0}$) ограниченного набора данных для C_P- и C_V-подходов перечислены в табл. 2 (Опт-3(C_P) и Опт-3(C_V) соответственно). Параметры, полученные при оптимизации Опт-3 по ограниченному набору данных, представлены в табл. S3 . Как видно из табл. 2, исключение из оптимизации определенных данных (PVT-данные при высоких температурах и давлениях) лишь незначительно ухудшает их расчет для обоих подходов с небольшим улучшением описания оставшихся данных. При этом качество оценки исключенных еще в оптимизации Опт-2 данных также меняется незначительно.

Вариант (в). При использовании C_P-варианта УС не удалось получить разумную оценку адиабатического модуля упругости, в то время как C_V-подход позволил рассчитать значения K_S при атмосферном давлении вплоть до температур 700–800 K.

На основании изложенного сделан вывод о том, что C_P-вариант обладает лучшими интерполяционными свойствами, а C_V-подход – экстраполяционными свойствами.

ЗАКЛЮЧЕНИЕ

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

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

Выявленные закономерности получены только на одном тестовом объекте – оксиде магния. О том, насколько они универсальны, можно будет судить при расширении круга изученных систем.