1. На главную
  2. Электронные версии
  3. Теплофизика высоких температур
  4. T. 58, № 1, 2020

Теплофизика высоких температур, 2020, T. 58, № 1, стр. 61-75

Компьютерное моделирование никеля и учет электронных вкладов в методе молекулярной динамики

Д. К. Белащенко *

Московский институт стали и сплавов
Москва, Россия

* E-mail: dkbel75@gmail.com

Поступила в редакцию 17.01.2019
После доработки 29.03.2019
Принята к публикации 16.05.2019

Полный текст (PDF)

Аннотация

Предложены два новых потенциала модели погруженного атома для никеля с учетом и без учета теплового электронного вклада в энергию. Параметры потенциалов найдены по свойствам никеля на изобаре p = 0 и в условиях ударного сжатия до давлений ~760 ГПа. Наилучшее согласие с экспериментом получено при учете электронного вклада в энергию. Расчетная адиабата Гюгонио хорошо согласуется с фактической. Учет электронного вклада сильно уменьшает температуру на адиабате и увеличивает “холодное” давление. Расчетная линия плавления моделей никеля плавно поднимается до 4518 К при 300 ГПа и почти нечувствительна к наличию электронного вклада. Определены координаты участка плавления на ударной адиабате (начало – 275.8 ГПа и 4422 К, конец – 297.6 ГПа и 4499 К). Предложена фазовая диаграмма никеля с областью устойчивой ОЦК-фазы при давлениях выше 110–130 ГПа. Приведены таблицы энергии и давления моделей при степенях сжатия до 1.8182.

ВВЕДЕНИЕ

Интерес к исследованию никеля обусловлен его важной ролью в металлургии, материаловедении и планетологии. Критические параметры никеля известны не точно. Согласно [1], критическая температура Tc лежит в интервале 5200–11 000 К, а критическая плотность dc в интервале от 0.60 до 2.3 г/см3. Молекулярно-динамические оценки [2] дают Tc = 9460 ± 20 К, dc = 2.560 ± ± 0.100 г/см3 и давление pc = 1.08 ± 0.01 ГПа. Согласно [3], стандартная теплота атомизации никеля равна $\Delta H_{{298}}^{0}$ = 431 кДж/моль. Отсюда энергия никеля в стандартных условиях по отношению к покоящимся удаленным атомам равна U298 = –[$\Delta H_{{298}}^{0}$ – (5/2)RT] = –424.8 кДж/моль. Энергию при повышенных температурах можно рассчитать по данным работ [4–8], где приведены значения энтальпии жидкого никеля при температурах до 3156 К. Так, при 1728 К получается U1728 = –359.6 кДж/моль. Данные по теплоемкости и энтальпии до 4300 К приведены в [9].

Плотность d жидкого никеля при температурах до 2030 К измерена в [10], до 2423 К – в [1] (d = (9.81–10.80) × 10–4T г/см3), до 2500 К в обзоре [11] и до 4258 К в [9] (d = (8.9817–6.4406) × × 10–4T г/см3, метод импульсного нагрева). В [9] измерена также скорость звука. Изотермический модуль всестороннего сжатия при 298 К равен KT = 180 ГПа [3] (и 168.2 ГПа в [12]), а данные при 1730 К (KT ≈ 99.5 ГПа [13]), видимо, завышены. Структура жидкого никеля исследована дифракционными методами при 1773–2023 К [14, 15] и 1873 К [16] (рис. 1). Структурные данные при 1773 К пересчитаны с помощью метода подавления ложных осцилляций парной корреляционной функции (ПКФ) [17]. Данные по ударному сжатию компактного никеля опубликованы в [18–20] и в обзоре [21].

Рис. 1.

Парные корреляционные функции (ПКФ) никеля: 1 – 1773 К, данные [14] с фурье-преобразованием по методу подавления ложных осцилляций [17]; 2 – 1873, [16]; 3 – 2023, [14]; маркеры – МД-данные, штриховые линии – дифракционные данные.

Опубликовано большое число работ по моделированию кристаллического никеля методом молекулярной динамики (МД). В ранних работах по моделированию жидкого никеля применялись потенциалы, составленные по известным дифракционным данным о структуре с помощью уравнения Перкуса–Йевика [22] или Борна–Грина–Боголюбова [23]. Этим методом можно достичь хорошего согласия ПКФ-модели и дифракционной ПКФ, однако обычно энергия модели отличается от фактической в несколько раз. Впоследствии были предложены потенциалы модели погруженного атома (Embedded Atom Model – EAM) для кристаллического [24–31] и жидкого [32–35] никеля. В [34] получено хорошее согласие структурного фактора модели жидкого никеля с дифракционным, но в целом потенциалы, предложенные для ГЦК-никеля, недостаточно хороши для жидкого металла.

В [26] применен метод ab initio для построения нескольких кристаллических структур никеля и для расчета параметров соответствующего потенциала ЕАМ. Полученный потенциал позволил правильно описать основные физические свойства кристаллического никеля, однако этот потенциал завышает плотность жидкого никеля при 2500–3600 К на ~1 г/см3 и занижает энергию на 10–20 кДж/моль. Кроме того, плотность жидкости в моделях оказывается выше, чем у ГЦК-никеля на ~0.2 г/см3.

В работах [33, 36] для моделирования жидкого никеля методом МД применен потенциал ЕАМ в форме [37]. Этот многочастичный потенциал имеет вид

(1)
$U = \sum\limits_i {\Phi ({{\rho }_{i}})} + \sum\limits_{i < j} {\varphi ({{r}_{{ij}}}).} $

Здесь r – расстояние в Å, Φ(ρi) – потенциал погружения i-го атома, зависящий от эффективной электронной плотности ρ в месте нахождения центра атома. Вторая сумма по парам атомов содержит обычный парный потенциал φ(r). Эффективная электронная плотность в точке нахождения атома создается окружающими атомами и определяется по формуле

(2)
${{\rho }_{i}} = \sum\limits_j {\psi ({{r}_{{ij}}})} ,$
где ψ(rij) – вклад в электронную плотность атома i от соседа номер j. В случае некристаллических металлов с нормальной или повышенной плотностью значения ρi у различных атомов близки к среднему и их разброс невелик (порядка нескольких процентов). Если приближенно считать все значения ρi одинаковыми, т.е. ρi = 〈ρ〉, то эффективная парная сила между атомами i и j равна
(3)
${{F}_{{ij}}}(r) = - 2{{\left. {\frac{{\partial {{{\Phi }}_{i}}}}{{\partial \rho }}\frac{{\partial \psi }}{{\partial r}}} \right|}_{{{{r}_{{ij}}}}}} - \frac{{\partial \varphi }}{{\partial r}},$
где Φi(ρ) – значение потенциала погружения на атомах. Если при каких-то условиях ∂Φ/∂ρ = 0, то межчастичные силы определяются только парным вкладом φ(r).

В [33, 36] рассчитаны свойства никеля на бинодали до 4273 К, а также в условиях ударного сжатия при давлениях до 764 ГПа, в двух вариантах: без учета тепловых вкладов коллективизированных электронов в энергию и давление (потенциал ЕАМ-1) и с учетом этих вкладов (потенциал ЕАМ-2). В [36] удалось рассчитать свойства никеля при высоких давлениях и температурах до 23 500 К, однако температуры плавления при высоких давлениях оказались заниженными.

В настоящей работе предложен новый потенциал ЕАМ для никеля, который позволяет улучшить результаты расчета свойств никеля и его линии плавления при высоких давлениях.

МЕТОДИКА ИССЛЕДОВАНИЯ

Парный вклад в потенциал ЕАМ. В [33, 36] парный вклад φ(r) в потенциал ЕАМ рассчитан с помощью алгоритма Шоммерса [38]. Полученный парный потенциал проходит через минимум около 3.90 Å, в то время как среднее межчастичное расстояние близко к 2.50 Å. При таких потенциалах кристаллическая решетка модели часто оказывается неустойчивой и разрушается при повышении температуры (как, например, в случае Nb и V с потенциалами [39]). Поэтому в настоящей работе парный вклад в потенциал ЕАМ аппроксимирован с использованием функции Морса:

(4)
$\begin{gathered} {{\varphi }_{{\text{M}}}}(r,\varepsilon ,\alpha ,{{r}_{0}}) = \varepsilon \{ {\text{exp[}}{\kern 1pt} --{\kern 1pt} 2\alpha ({r \mathord{\left/ {\vphantom {r {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}}--1) - \\ --\,\,2{\text{exp[}}{\kern 1pt} --{\kern 1pt} \alpha ({r \mathord{\left/ {\vphantom {r {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}}--1){\text{]}}\} , \\ \end{gathered} $
а при r ≤ rc

(5)
$\begin{gathered} \varphi (r) = {{\varphi }_{{\text{M}}}}(r,{{\varepsilon }_{1}},{{\alpha }_{1}},{{r}_{0}}) + {{\varepsilon }_{2}}{\text{exp[}}{\kern 1pt} --{\kern 1pt} {{\alpha }_{2}}{{({{r}_{2}}--r)}^{2}}{\text{]}} + \\ + \,\,{{\varepsilon }_{3}}{{({{r}_{c}}--r)}^{2}}--{{\varphi }_{{\text{M}}}}({{r}_{c}},{{\varepsilon }_{1}},{{\alpha }_{1}},{{r}_{0}}). \\ \end{gathered} $

Ранее потенциал Морса (4) для описания парного вклада никеля выбирался в [2] (с параметрами из [40]: ε = 0.4205 эВ, α = 3.947 и r0 = 2.78 Å). Дополнительные слагаемые в (5) позволяют получить правильное значение координаты не только первого, но и второго пиков ПКФ, а также близкое к нулю давление модели. Значения параметров ε1, ε2, ε3, α1, α2, r0, r2, rc приведены в табл. 1. Этот парный вклад в потенциал ЕАМ-1 никеля (с минимумом около 2.63 Å) показан на рис. 2. У жидкости с таким потенциалом (не ЕАМ!) и с реальной плотностью никеля давление при 1773 К близко к нулю и ПКФ практически совпадает с ПКФ реального никеля (рис. 1). Стандартное отклонение между этими ПКФ (“невязка”) при 1773 К составляет всего Rg = 0.042.

Таблица 1.  

Параметры потенциала ЕАМ никеля

Параметр Значение, эВ, Å Параметр Значение Параметр Значение, эВ Параметр Значение, эВ Параметр Значение, эВ
Парный вклад Потенциал погружения
ε1 0.150 ρ1 0.89 a1 –2.594400 b1 – c1 1.22920
ε2 0.050 ρ2 0.78 a2 –2.579527 b2 –0.270424 c2 –3.520
ε3 –0.00189 ρ3 0.69 a3 –2.592372 b3 0.503976 c3 8.000
α1 4.20 ρ4 0.62 a4 –2.572930 b4 –0.936024 c4 –2.100
α2 7.78 ρ5 0.46 a5 –2.517698 b5 –0.642024 c5 1.000
r0 2.75 ρ6 1.15 a6 –2.389374 b6 –0.962024 c6 1.000
r2 3.00 ρ7 1.935 a7 –2.566743 b7 0.368760 c7 1.340
r3 5.50 M 1.90 a8 –1.431292 b8 2.416341 c8 (ЕАМ–A) 2.710
rc 7.05 N 2.00 – – – – c8 (ЕАМ-B) 3.480
p1 0.001559 – – – – – – – –
Рис. 2.

Парный вклад в потенциал ЕАМ-1 никеля.

Эффективная электронная плотность. Эффективная электронная плотность ЕАМ выбрана в форме

$\psi (r) = {{p}_{1}}{{({{r}_{3}}--r)}^{{3.6}}}{{H}_{1}}({{r}_{3}}--r),$
где H1(z) = 1 при z ≥ 0 и H1(z) = 0 при z < 0. Такая форма обеспечивает достаточно высокую чувствительность средней эффективной плотности электронов на атомах 〈ρ〉 к изменениям плотности никеля при нагревании. Значения p1 и r3 приведены в табл. 1. При таких значениях параметров эффективная электронная плотность ρi = $\sum\nolimits_j {\psi ({{r}_{{ij}}})} $ на атомах жидкого никеля при 1773 К и давлении p ≈ 0 близка к единице (1.001 ± 0.072).

Потенциал погружения. Потенциал погружения выбран в следующей кусочно-непрерывной форме:

$\begin{gathered} \Phi \left( \rho \right) = {\text{[}}{{a}_{1}} + {{c}_{1}}{{\left( {\rho --1} \right)}^{2}}{\text{]}}{{H}_{{12}}}\left( {\rho ,{{\rho }_{1}},{{\rho }_{6}}} \right) + \\ + \,\,[{{a}_{2}} + {{b}_{2}}\left( {\rho --{{\rho }_{1}}} \right) + {{c}_{2}}{{\left( {\rho --{{\rho }_{1}}} \right)}^{2}}]{{H}_{{12}}}\left( {\rho ,{{\rho }_{2}},{{\rho }_{1}}} \right) + {\text{ }} \\ {\text{ + }}\,\,[{{a}_{3}} + {{b}_{3}}\left( {\rho --{{\rho }_{2}}} \right) + {{c}_{3}}{{\left( {\rho --{{\rho }_{2}}} \right)}^{2}}]{{H}_{{12}}}\left( {\rho ,{{\rho }_{3}},{{\rho }_{2}}} \right) + \\ + \,\,[{{a}_{4}} + {{b}_{4}}\left( {\rho --{{\rho }_{3}}} \right) + {{c}_{4}}{{\left( {\rho --{{\rho }_{3}}} \right)}^{2}}]{{H}_{{12}}}\left( {\rho ,{{\rho }_{4}},{{\rho }_{3}}} \right) + \\ + \,\,[{{a}_{5}} + {{b}_{5}}\left( {\rho --{{\rho }_{4}}} \right) + {{c}_{5}}{{\left( {\rho --{{\rho }_{4}}} \right)}^{2}}]{{H}_{{12}}}\left( {\rho ,{{\rho }_{5}},{{\rho }_{4}}} \right) + \\ + \,\,[{{a}_{6}} + {{b}_{6}}\left( {\rho --{{\rho }_{5}}} \right) + {{c}_{6}}{{\left( {\rho --{{\rho }_{5}}} \right)}^{2}}] \times \\ \times \,\,[{{2\rho } \mathord{\left/ {\vphantom {{2\rho } {{{\rho }_{5}}}}} \right. \kern-0em} {{{\rho }_{5}}}}--{{\left( {{\rho \mathord{\left/ {\vphantom {\rho {{{\rho }_{5}}}}} \right. \kern-0em} {{{\rho }_{5}}}}} \right)}^{2}}]{{H}_{{12}}}\left( {\rho ,0,{{\rho }_{5}}} \right) + \\ + \,\,[{{a}_{7}} + {{b}_{7}}\left( {\rho --{{\rho }_{6}}} \right) + {{c}_{7}}{{\left( {\rho --{{\rho }_{6}}} \right)}^{m}}]{{H}_{{12}}}\left( {\rho ,{{\rho }_{6}},{{\rho }_{7}}} \right) + \\ + \,\,[{{a}_{8}} + {{b}_{8}}\left( {\rho --{{\rho }_{7}}} \right) + {{c}_{8}}{{\left( {\rho --{{\rho }_{7}}} \right)}^{n}}]{{H}_{{12}}}\left( {\rho ,{{\rho }_{7}},\infty } \right). \\ \end{gathered} $

Здесь H12(z, a, b) = 1 при a ≤ z ≤ b и H12 = 0 в остальных случаях. Очень полезная особенность такого выбора функции заключается в том, что производная dΦ(ρ)/dρ обращается в нуль при ρ = 1, так что при условии 〈ρ〉 = 1 движением атомов управляет только парный потенциал φ(r) согласно (3).

Параметры a1 и c1–c6 подобраны по значениям энергии жидкого никеля при 1773 К и по зависимости плотности жидкого никеля от температуры [11]. В частности, параметр a1 выбран таким образом, чтобы вблизи точки плавления энергия жидкого никеля с учетом электронного вклада равнялась экспериментальной величине –359.6 кДж/моль. Параметры ai и bi (i = 1–8) найдены из условий непрерывности функции Φ(ρ) и ее первой производной. Значения параметров приведены в табл. 1.

Для расчета параметров c7 и c8 потенциала погружения использовались результаты ударного сжатия компактного никеля, опубликованные в [18–21]. Ударная адиабата, построенная по этим данным, показана на рис. 3. Она хорошо аппроксимируется выражением

$\begin{gathered} p = 1233.471{{Z}^{3}}--3923.964{{Z}^{2}} + \\ + \,\,4439.137Z--1752.968, \\ \end{gathered} $
где V – мольный объем, Z = V0/V и V0 = = 6.6133 см3/моль. Энергия на адиабате U рассчитывается по уравнению [41]
$U--{{U}_{0}} = 0.5\left( {p + {{p}_{0}}} \right)\left( {{{V}_{0}}--V} \right),$
где p0 и U0 – исходные давление и энергия. Параметры c7 и c8 подобраны так, чтобы получить совпадающие с фактическими значения давления и энергии МД-модели на адиабате Гюгонио, причем либо без учета тепловой энергии электронов (потенциал ЕАМ-A), либо с учетом (ЕАМ-B). Методика расчета описана, например, в [33, 36]. Найденные подгонкой значения параметров потенциала погружения приведены в табл. 1. Потенциал погружения показан на рис. 4 для обоих случаев учета/неучета электронных вкладов в энергию никеля. Эти потенциалы отличаются только значением параметра c8 и приводят к различным результатам лишь при ρ > ρ7 = 1.935.

Рис. 3.

Адиабата Гюгонио никеля: 1 – эксперимент [18–20]; 2 – МД-расчет с потенциалом ЕАМ-A без учета электронных вкладов; 3 – ЕАМ-B, с учетом электронных вкладов в энергию; V0 = 6.6133 см3/моль.

Рис. 4.

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

Учет тепловых свойств электронов никеля. Электронный вклад в теплоемкость никеля рассчитан в [42–44] (в единицах Дж/(м3 К)) по данным для плотности d-состояний. Эти данные пересчитаны на размерность Дж/(моль К) в предположении, что плотность никеля постоянна и равна 8.875 г/см3 (мольная масса 58.693). По данным [42, 43] в интервале 0 < T < 105 К электронная теплоемкость аппроксимируется выражением

(6)
$\begin{gathered} {{C}_{{eT}}} = --8.986638 \times {{10}^{{--12}}}{{T}^{3}} + 8.538552 \times {{10}^{{--8}}}{{T}^{2}} + \\ + \,\,1.345664 \times {{10}^{{--3}}}T + 2.027078, \\ \end{gathered} $
а в интервале 0–5 × 104 К выражением

(7)
$\begin{gathered} {{C}_{{eT}}} = 1.146241 \times {{10}^{{--21}}}{{T}^{5}}--1.656827 \times {{10}^{{--16}}}{{T}^{4}} + \\ + \,\,8.530214 \times {{10}^{{--12}}}{{T}^{3}}--1.957821 \times {{10}^{{--7}}}{{T}^{2}} + \\ + \,\,2.690477 \times {{10}^{{--3}}}T + 0.5583407. \\ \end{gathered} $

Тепловая энергия электронов EeT рассчитывается интегрированием этих выражений. Ее график показан на рис. 5.

Рис. 5.

Тепловая энергия электронов никеля по данным [42, 43].

Для расчетов существенен способ учета теплового давления электронов. Давление электронов peT рассчитывается по уравнению peTV = λEeT (V – объем металла). Поведение электронного коэффициента Грюнайзена λ обсуждалось в литературе неоднократно. В модели свободных электронов λ = 2/3 [45]. Для невзаимодействующих электронов – 0 < λ < 2/3. Случаи более сложного поведения электронов рассмотрены в [46]. Расчеты электронного давления для Al, Au, Cu, K, Ni, Ta, Ti, W проведены, например, в [44, 47, 48], в частности для случая сверхскоростного лазерного нагрева металла. Для использования в МД-моделировании требуется подходящая параметризация этих данных.

Метод моделирования. Моделирование проводилось методом МД с алгоритмом Л. Верле. Модели содержали 2000 (ОЦК) или 2048 (ГЦК) частиц в основном кубе. Свойства моделей определялись в режиме NVT, а расчеты линии плавления – в режиме NpT при ступенчатом нагреве. Длина изотермических прогонов составляла обычно 40 000–50 000 шагов по времени, длина шага Δt = 0.01t0, где внутренняя единица времени t0 = 7.800 × 10–14 с. Температура плавления находилась методом отогрева.

РЕЗУЛЬТАТЫ РАСЧЕТОВ

В настоящей работе рассмотрены два варианта поведения электронов. Первый – в предположении, что электроны локализованы на атомах (как, например, у урана [36]) и не дают вклада ни в энергию, ни в давление (потенциал ЕАМ-A), а второй – в предположении, что электроны дают вклад в энергию (согласно [42, 43]), но не в давление (т.е. λ = 0, потенциал ЕАМ-B). Сравнение этих вариантов позволяет оценить роль тепловых свойств электронов при применении метода МД.

Свойства моделей никеля на бинодали. Результаты расчетов с учетом тепловой энергии электронов приведены в табл. 2. Электронная энергия EeT (формулы (6), (7), рис. 5) добавлена к значениям энергии моделей UMD. При давлениях меньше ~200 ГПа потенциалы ЕАМ-A и ЕАМ-B дают совпадающие результаты. Как видно из табл. 2, потенциал ЕАМ, подобранный для жидкого никеля, неплохо описывает энергию ГЦК-никеля при 298 К (отличие всего на 4.4 кДж/моль), но занижает плотность на 1.6%. Энергия ГЦК-никеля в этих условиях ниже, чем у ОЦК-никеля, так что структура ГЦК более устойчива, чем ОЦК. Теоретические расчеты [26] дают близкие результаты, но завышают разницу объемов и энергий ГЦК- и ОЦК-никеля при 0 К. Расчетная плотность жидкости очень хорошо согласуется с экспериментом вдоль бинодали до 2573 К, а при более высоких температурах согласие ухудшается, возможно, из-за меньшей точности данных [9], полученных методом импульсного нагрева (они явно сдвинуты вверх примерно на 0.45 г/см3 по сравнению с данными обзора [11]). Хорошее согласие с экспериментом по энергии жидкости получается только при учете тепловой энергии электронов. Структурные характеристики моделей хорошо согласуются с дифракционными данными (рис. 1), и невязки Rg при 1773–2023 К равны 0.040–0.045.

Таблица 2.  

Свойства МД-моделей никеля с потенциалом ЕАМ-B при p ≈ 0

T, К d, г/см3 〈ρ〉b Rg U, кДж/моль D × 105, см2/с KT, ГПа Теплоемкость, Дж/(моль К) u, м/с
МД [1, 11, 41] UMD EeT [42, 43] UMD + EeT Uexp
[3, 6, 8, 9] 		CeT $C_{V}^{{\text{d}}}$ $C_{p}^{{\text{d}}}$
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
ОЦК
298 8.739 – 1.1378 – –417.99 0 –417.99 –424.8с – 191 0 25.6 25.8 4697
ГЦК
298 8.736 8.875 1.1374 – –420.38 0 –420.38 –424.8c – 198 0 28.8 29.0 4798
1100 8.541 8.59 1.1088 – –398.79 2.42 –396.37 –399.8 – 169.7 3.56 25.6 27.3 4583
Жидкий никель
1730 7.848 7.858 1.0064 – –364.46 5.01 –359.45 –359.6 6.15 77.69 4.50 28.50 34.18 3446
1773 7.812 7.816 1.0015 0.045 –362.98 5.20 –357.78 –357.9 6.78 79.99 4.49 29.79 35.49 3582
1873 7.748 7.716 0.9895 0.042 –359.89 5.68 –354.21 –354.0 7.21 62.54 4.66 30.18 37.08 3149
2023 7.619 7.568 0.9712 0.040 –354.86 6.41 –348.45 –348.2 9.13 – – – – –
2073 7.582 7.518 0.9655 – –353.31 5.67 –347.64 –346.3 10.1 36.94 5.00 27.40 39.82 2661
2273 7.340 7.320 0.9318 – –345.90 7.72 –338.18 –338.6 12.8 – – – – –
2573 6.962 7.025 0.8710 – –335.10 9.43 –325.67 –327.1 16.6 26.18 5.85 27.33 42.73 2425
3073 6.546 7.002a 0.8326 – –320.10 12.60 –307.50 –307.8 25.4 29.00 6.72 30.84 33.48 2193
3573 6.262 6.680a 0.7727 – –306.48 16.17 –290.31 –286.6а 31.4 – – – – –
4073 6.003 6.358a 0.7894 – –293.21 20.15 –273.06 –265.4а 38.2 – – – – –
4273 5.812 6.230a 0.7077 – –286.22 21.85 –264.37 –256.9а 45.6 – – – – –

a Данные [9] при давлении 0.2 ГПа, остальные при p ≅ 0. b Стандартное отклонение между значениями ρ на атомах модели растет сверху вниз от 0.025 до 0.093. с Данные [3]. d С учетом СeT по (6) и (7).

В табл. 2 приведены также значения теплоемкостей CV и Cp, модуля всестороннего сжатия KT и скорости звука u, рассчитанные по МД-данным при нескольких температурах. Экспериментальное значение теплоемкости Cp в интервале 1730–2000 К равно 38.9 Дж/(моль К) [6] и близко к данным табл. 2 с учетом электронного вклада. Модуль всестороннего сжатия KT при 298 К (198 ГПа) близок к фактическому (186.7 [49], 180 ГПа [3]) и к значениям, полученным в других работах по моделированию никеля (178.7 [50], 196.1 [51]). Модуль KT и скорость звука u модели ГЦК-никеля при 298 К неплохо согласуются с экспериментом (KT = 180 ГПа и u = 4970 м/с [3]). Экспериментальные значения KT жидкого никеля, найденные по данным [9] для скорости звука, выше приведенных в табл. 2 примерно на 20%. Эти расхождения коррелируют с завышением плотности никеля в [9] на ~0.45 г/см3. Кроме того, расчетная скорость звука в табл. 2 убывает с температурой гораздо быстрее, чем в [9].

Коэффициенты самодиффузии D, рассчитанные по зависимости от времени среднего квадрата смещений частиц (табл. 2), аппроксимируются уравнением

$D = 5.344 \times {{10}^{{--12}}}{{T}^{{2.1887}}}$

с разбросом около 10%. Экспериментальные данные по самодиффузии никеля [52, 53] (метод неупругого рассеяния нейтронов) дают 3.8 × 10–5 при 1795 К. Авторы [54] получили методом ab initio значение D = 4.4 × 10–5 при 1850 К, что несколько меньше настоящих данных. В [32] (метод МД) при 1773 К также получено более низкое значение D = 4.8 × 10–5. Однако МД-моделирование в [55] (с потенциалом Саттона–Чена) дало значения коэффициента самодиффузии при обычном давлении от 6.69 × 10–5 при 1773 К до 11.7 × 10–5 при 2273 К. Эти цифры очень близки к приведенным в табл. 2.

Расчет адиабаты Гюгонио моделей никеля. Свойства моделей с потенциалом ЕАМ-B при условиях ударного сжатия приведены в табл. 3. На рис. 3 показаны расчетные точки ударной адиабаты, полученные с применением потенциалов ЕАМ-A и ЕАМ-B. Согласие расчетных данных с фактическими здесь очень хорошее за исключением области, где p < 80 ГПа, поскольку потенциалы ЕАМ-A и ЕАМ-B, видимо, недостаточно точны для описания кристаллических фаз. Согласие фактических давлений на адиабате и расчетных давлений при Z = V0/V > 1.20 хорошее; фактические энергии и энергии моделей также довольно близки (табл. 3). Разница в температурах THug вдоль адиабат для обоих потенциалов (ЕАМ-A и ЕАМ-B) появляется при давлениях на адиабате свыше 100 ГПа, но в интервале 100 < p < 300 ГПа она невелика и составляет 200–300 К. При более высоких давлениях разница расчетных THug быстро увеличивается.

Таблица 3.  

Свойства моделей никеля на ударной адиабате. Потенциал ЕАМ-B

Z d, г/см3 p, ГПа U2 – U1,
кДж/моль		 ТHug, К 〈ρ〉 d U, кДж/моль pMD, ГПа, модель
ЕАМ-A1  ЕАМ-B1 [56]1 [57]e UMD EeT
(6), (7) 		UMD + EeT U298 + U2 – U1
1.00a 8.736c 0 0.00 298 298 298 298 1.137 –420.38 0 –420.38 –424.8 0.0
1.111a 9.861 21.8 7.2 300 300 310 – 1.331 –410.3 0 –410.3 –417.6 32.4
1.176a 10.441 48.8 24.2 310 310 550 – 1.432 –396.5 0 –396.5 –400.6 53.2
1.250a 11.094 78.1 51.7 350 350 1040 905 1.546 –373.7 0.12 –373.6 –373.1 80.0
1.333a 11.833 114.7 94.8 820 800 1740 1385 1.677 –331.0 1.40 –329.6 –330.0 116.7
1.429a 12.678 172.9 171.6 2200 2040 3130 2340 1.831 –255.6 6.48 –249.1 –253.2 171.6
1.429b 12.678 172.9 171.6 2300 2040 3130 – 1.829 –259.9 6.48 –253.4 –253.2 170.8
1.538 13.654 277.0 320.5 4600 4450 6170 4390 2.007 –126.1 23.39 –102.7 –104.3 274.6
1.667 14.792 460.2 608.7 12 250 8560 12 520 7800 2.217 104.0 72.04 176.0 183.9 465.9
1.818 16.136 765.7 1139.4 30 300 17 820 – – 2.534 485.4 230.0 715.4 714.6 763.2

a ГЦК-Ni. b ОЦК-Ni. c При p = 0. d Стандартное отклонение растет сверху вниз от 0.015 до 0.11. e Интерполяция.

При учете электронных вкладов (ЕАМ-B) температура на адиабате заметно ниже, чем рассчитанная аналитически в [56]. При давлениях ниже 277 ГПа расчеты [57] завышают температуру на адиабате, а при более высоких давлениях занижают ее по отношению к результатам с ЕАМ-B (табл. 3).

Холодное давление никеля. В работах [58, 59] проведено статическое сжатие ряда металлов в алмазной ячейке при 298 К, в том числе никеля до давления 156 ГПа. Используя потенциалы ЕАМ-A и ЕАМ-B, можно рассчитать давление моделей никеля при температурах 0 или 298 К. Давления при этих температурах близки друг к другу и даже при максимальной степени сжатия отличаются всего на 1.5 ГПа. В табл. 4 приведены значения “холодного давления” моделей ОЦК- и ГЦК-никеля при 298 К в зависимости от плотности. При p < < 200 ГПа эти данные в целом хорошо согласуются с результатами статического сжатия никеля [58] и с данными для холодного давления, полученными обработкой экспериментов по ударному сжатию [56, 60], а при p > 200 ГПа результаты расчетов с потенциалами ЕАМ существенно выше данных [60]. При p > 200 ГПа учет электронных вкладов приводит к увеличению холодного давления при заданной плотности (на 9 ГПа при 14 г/см3 и на 56 ГПа при 16 г/см3). В [36] принято, что CeT = $C_{{eT}}^{0}$Z –2/3 (как в модели свободных электронов), что привело к понижению холодного давления (до 459 ГПа при Z = 1.8182). При давлениях ~150 ГПа квантово-механические расчеты [58] дают разброс в ~30 ГПа по отношению к статическим данным, а расчеты [50] с потенциалом ЕАМ (Саттона–Чена) завышают его на ~40 ГПа (в [58] проведены как данные измерений сжимаемости никеля при 300 К, так и результаты квантово-механических расчетов сжимаемости).

Таблица 4.  

Давление никеля и моделей ОЦК- и ГЦК-никеля при 298 К без учета (1) и с учетом (2) тепловых вкладов электронов

Система d, г/см3
8.8 9.0 9.5 10.0 11.0 12.0 13.0 14.0 15.0 16.0
p, ГПа
ОЦК-1 3.1 9.2 23.2 36.8 68.8 115.7 185.3 282.2 408.1 561.8
ОЦК-2 2.9 10.2 25.5 38.7 67.6 112.4 184.6 291.4 436.0 617.2
ГЦК-1 0.7 4.9 18.7 35.9 75.9 121.7 181.0 272.1 – –
ГЦК-2 0.5 3.6 16.8 35.0 76.9 121.8 180.9 283.7 – –
Cтатическое сжатие, [58] 1.4 4.9 15.9 30.1 68.2 119.0 182.9 – – –
Интерполяция [60] – 4.4 22.7 38.8 72.1 117.0 181.1 267.6 374.6 495.4
[61] – 5.5 15.5 28.3 62.6 109.4 169.7 244.3 334.2 440.3
Квантово-механический расчет, [58] – – 4.0 17.2 53.5 103.4 166.8  244a – –
[50] – 0.6 10.5 25.2 67.0 125.7 207.2 324.2a – –

a Экстраполяция.

Относительную устойчивость решеток ОЦК и ГЦК можно оценить по величинам их энергий Гиббса G = U + pV – TS (S – энтропия). При 298 К и давлении ~50 ГПа теплоемкости Cv ОЦК- и ГЦК-фаз практически совпадают (соответственно 25.11 и 25.13 Дж/(моль К)), так что разницу энтропий этих фаз можно считать несущественной. По данным табл. 5 получается, что при 298 К структура ГЦК устойчивей, чем ОЦК, примерно до 52 ГПа (ее энергия Гиббса ниже), а при более высоких давлениях устойчива структура ОЦК. Переход от ГЦК-структуры никеля, как наиболее устойчивой, к ОЦК-структуре при сжатии металла виден также по результатам расчетов в [26], проведенным как методом ab initio, так и на основе потенциала ЕАМ. Следовательно, на фазовой диаграмме никеля должна наблюдаться область устойчивой ОЦК-фазы, с границей ГЦК–ОЦК при 298 К около 52 ГПа. Вторую точку на этой границе можно получить при анализе линии плавления моделей никеля (см. ниже).

Таблица 5.  

Плотность и энергия моделей ОЦК- и ГЦК-никеля при 298 К без учета (1) и с учетом (2) тепловых вкладов электронов

Система p, ГПа
0 50 100 150 200 250 300 400 500 600
d, г/см3
ОЦК-1 8.713 10.564 11.655 12.512 13.238 13.728 14.159 14.94 15.62 16.24
ОЦК-2 8.713 10.564 11.655 12.512 13.236 13.678 14.060 14.75 15.36 15.92
ГЦК-1 8.730 10.354 11.533 12.494 13.261 13.773 14.222 – – –
ГЦК-2 8.730 10.354 11.533 12.493 13.252 13.716 14.108 – – –
  Энергия, кДж/моль
ОЦК-1 –417.4 –392.9 –356.0 –313.9 –269.7 –234.5 –199.0 –124.3 –47.3 30.6
ОЦК-2 –417.4 –392.9 –356.0 –313.9 –269.8 –238.1 –206.2 –138.4 –67.4 –5.4
ГЦК-1 –420.2 –398.7 –356.9 –309.1 –262.5 –226.1 –189.9 – – –
ГЦК-2 –420.2 –398.7 –356.9 –309.1 –263.1 –230.0 –197.7 – – –

Зависимость холодного давления от степени сжатия при 298 К при использовании потенциала ЕАМ-B описывается уравнением

$\begin{gathered} {{p}_{{298}}} = --1515.0{{Z}^{4}} + 9526.2{{Z}^{3}}--20{\kern 1pt} 927{{Z}^{2}} + \\ + \,\,19{\kern 1pt} 810Z--6895.9. \\ \end{gathered} $

На рис. 6 показаны графики холодного давления, полученные статическим методом [58], обработкой результатов ударного сжатия [56, 60] (линии 1 и 2) и рассчитанные с использованием потенциалов ЕАМ-A и ЕАМ-B. Они заметно различаются лишь при Z > 1.5. Учет электронных вкладов приводит к росту расчетного холодного давления (линия 4).

Рис. 6.

“Холодное давление” никеля: 1 – данные [60] при 0 К; 2 – [56], 0 К; 3 – авторские данные при 298 К, потенциал ЕАМ-A; 4 – при 298 К, ЕАМ-B; 5 – статическое сжатие [58].

Линия плавления моделей никеля. Температура плавления Tm определялась методом отогрева при ступенчатом нагревании модели [36]. При каждой температуре проводилась серия изотермических прогонов длиной по 50 000 шагов по времени при наблюдении за средним квадратом смещений частиц и за максимальной величиной структурного фактора:

$S(q) = \frac{1}{N}{{\left| {\sum\limits_i {{{e}^{{ - iq{{R}_{i}}}}}} } \right|}^{2}}.$

Здесь N – число частиц модели, Ri – координаты атомов, q – вектор рассеяния излучения, осреднения по всем направлениям не проводится. При размере модели в несколько тысяч атомов максимальная величина S(q) кристаллической модели (в направлении векторов обратной решетки) близка к N, а после плавления убывает скачком до 20–30. С учетом флуктуаций температуры модели точку плавления можно детектировать с ошибкой ±5 К. Метод отогрева дает хорошие результаты при расчетах Tm для Li, Na, K, Rb, Cu, Ga, Pb, Fe [36] (согласие с экспериментом в основном в пределах 20–30 К). При этом он в техническом отношении гораздо удобнее метода моделирования двухфазной системы, который не всегда точен и приводит в ряде случаев к значительным ошибкам (например, дает отклонение в 90 К для температуры плавления никеля [62]). С другой стороны, метод детектирования плавления в процессе непрерывного нагревания модели кристалла приводит к значительному перегреву и завышению Tm [63]. Можно утверждать, что расхождения с экспериментом по температуре плавления особенно при высоких давлениях обусловлены не только спецификой метода моделирования, но и неполной адекватностью потенциала, который, как правило, недостаточно точен сразу для обеих конкурирующих фаз. Если теплота плавления в модели и скачок объема при плавлении заметно отличаются от фактических, то зависимость Tm от давления неизбежно получится неправильной. Вообще судить о качестве расчета температуры плавления следует не по одной точке, а по всей линии плавления Tm(p). В частности, хорошие результаты получены методом отогрева для линии плавления железа до давления 420 ГПа [64].

В табл. 6 приведены найденные методом отогрева температуры плавления моделей никеля при различных давлениях с учетом и без учета электронных вкладов. Температура плавления модели ГЦК-фазы при нулевом давлении равна 1656 К и немного меньше фактической Tm = 1728 К. Плотность в ГЦК-модели при 1660 К равна 8.267 г/см3, а плотность жидкой фазы 7.884, так что скачок объема при плавлении составляет ΔV = = 0.329 см3/моль, или 4.63%. Эта цифра хорошо согласуется с фактической величиной 4.70 ± ± 0.15% [10]. Теплота плавления по модели ΔH = 12.97 кДж/моль. Отсюда наклон линии плавления в этих условиях составляет dTm/dp = = TΔV/ΔH = 42.1 К/ГПа. Фактический начальный наклон линии плавления никеля равен 33 К/ГПа [65]. При давлениях 0–100 ГПа на линии плавления моделей устойчива ГЦК-фаза (температура плавления выше, чем у ОЦК), а при более высоких давлениях более устойчива ОЦК-фаза никеля. При ~126 ГПа температуры плавления ГЦК- и ОЦК-моделей становятся равными, определяя тем самым тройную точку при ~3500 К на линии плавления. Эти данные позволяют построить фазовую диаграмму никеля с дополнительной областью устойчивости ОЦК-фазы при высоких давлениях (рис. 7). Граница областей устойчивости ГЦК- и ОЦК-фаз должна идти от точки при 298 К и ~52 ГПа к тройной точке при 3500 К и ~125 ГПа. Авторы [66] также обнаружили признаки фазового перехода ГЦК–ОЦК при МД-моделировании ударного сжатия никеля.

Таблица 6.  

Зависимость температуры плавления моделей никеля от давления

Структура Метод p, ГПа
0 25 50 100 120 150 200 250 300 400 500 600
температура плавления, К
ОЦК ЕАМ-A 1437 – 2544 3231 3450 3706 4056 4356 4518 4843 5118 5332
ЕАМ-B 1437 – 2544 3231 3450 3706 4090 4332 4518 4806 5040 5281
ГЦК ЕАМ-A 1656 2332 2733 3283 – 3668 4018 4281 4456 – – –
ЕАМ-B 1656 2332 2733 3283 – 3668 4031 4318 4456 – – –
Рис. 7.

Фазовая диаграмма никеля (проект): γ – ГЦК-фаза, β – ОЦК-фаза, L – жидкость.

Очень важно то, что Tm данной фазы довольно слабо зависит от учета или неучета электронных вкладов (не только при давлениях до 300 ГПа, когда от этого мало зависит температура на адиабате, но и при давлениях до 600 ГПа, когда температура в двух вариантах может отличаться на несколько тысяч градусов). Различие значений Tm не превышает в рассматриваемых случаях 80 К при 500 ГПа (табл. 6). Кроме того, по данным автора, Tm алюминия при давлении 200 ГПа, рассчитанные без учета и с учетом тепловых вкладов электронов, равны 6792 и 6695 К (при разбросе статических данных разных авторов ~100 К), а в случае железа при давлении 400 ГПа они равны соответственно 7740 и 7300 К [64, 36 ], т.е. разница составляет от 1.4 до 5.7%. Следовательно, для МД-расчета линий плавления можно использовать потенциал ЕАМ, построенный по данным ударного сжатия без учета электронных вкладов – по крайней мере, в тех случаях, когда электронным давлением можно пренебречь. Этот метод фактически используется в большинстве МД-расчетов линии плавления металлов. Однако для систем с высоким электронным давлением требуется дополнительное исследование роли электронных вкладов.

На рис. 8 показаны результаты исследования температуры плавления никеля различными методами. В работах [29, 55, 64, 62, 71, 72 ] и ряде других выполнены теоретические расчеты методом ab initio и с использованием различных потенциалов ЕАМ. Полученные в них значения Tm выше рассчитанных здесь. Наоборот, в экспериментальных работах [67–69] (метод сжатия на алмазной наковальне – DAC) значения Tm ниже полученных здесь с потенциалами EAM-A и EAM-B, и только в [70] выше (на ~500 К при 125 ГПа). Наибольшие значения Tm получены с применением аналитической программы PANDA [71]. Причины расхождения экспериментальных данных подробно рассмотрены в [70].

Рис. 8.

Температура плавления никеля: 1 – настоящая работа, 2 – данные [55], 3 – [29], 4 – [62], 5 – [65, 67], 6 – [68, 69], 7 – [70], 8 – [57].

В [36] рассчитаны линии плавления никеля с предложенными там потенциалами ЕАМ-1 и ЕАМ-2. Найденная в [36] линия плавления сильно отличается от полученной выше и идет очень полого, занижая Tm. Причиной этого является, видимо, неудачная форма парного вклада в потенциал, найденного алгоритмом Шоммерса (см. выше). Добавление к потенциалу (4) функции Δ = 45.0(1.85 – r)2H1(1.85 – r) увеличивает крутизну отталкивательной ветви парного вклада в потенциал и увеличивает расчетную температуру плавления на 70–800 К.

Пересечение линии плавления с ударной адиабатой. Для более детального описания эффектов плавления при ударном сжатии зависимости свойств твердой и жидкой фаз в точках плавления (плотность, давление, температура) от степени сжатия интерполированы полиномами. Затем рассчитаны две ударные адиабаты: для твердой и жидкой фаз по отдельности. Метод расчета описан, например, в [33, 36]. Для каждой из этих адиабат определены точки пересечения с линией плавления.

На рис. 9 изображены кривые температуры на ветвях адиабаты Гюгонио (при учете электронного вклада) для твердой и жидкой фаз совместно с линией плавления. Точки пересечения линий 1–3 и 2–3 указывают на начало и конец плавления ОЦК-никеля при ударном сжатии. В данном случае плавление начинается при давлении 275.8 ГПа и температуре 4422 К, а заканчивается при давлении 297.6 ГПа и температуре 4499 К. Интервал плавления невелик, что и затрудняет его обнаружение в экспериментах по ударному сжатию. Положение точки плавления на ударной адиабате существенно зависит от способа расчета. Так, в [71] расчет с помощью аналитической программы PANDA дает начало плавления никеля при давлении 278 ГПа и окончание при 377 ГПа. В [57] ударная адиабата никеля рассчитана аналитически, и температура на ней выше приведенной в табл. 3. Плавление начиналось при давлении 350 ГПа (5900 К) и заканчивалось при 420 ГПа.

Рис. 9.

Определение температуры начала и конца плавления на ударной адиабате моделей никеля с потенциалом ЕАМ-B: 1 – для твердой фазы, 2 – для жидкой фазы, 3 – Tm.

Расчет свойств моделей никеля при высоких давлениях и температурах. С использованием потенциала ЕАМ-B построены модели и рассчитаны свойства никеля при степенях сжатия до Z = = 1.818 и температурах до 25 000 К. К значениям энергии построенных моделей добавлена тепловая энергия электронов по (6), (7). Расчетные энергии приведены в табл. 7, а давления – в табл. 8. Они согласуются с найденными в [36] при параметрах ударной адиабаты. Эти данные достаточны для расчета по обычным формулам остальных термодинамических свойств: теплоемкостей CV и Cp, производной (∂p/∂T)V, KT, u, коэффициента теплового расширения α и коэффициента Грюнайзена γ = (V/CV)(∂p/∂T)V. Расчеты выполнены для твердой и жидкой фаз отдельно с аппроксимацией полиномами различной степени. В данной работе приведены результаты только для модуля всестороннего сжатия и коэффициента Грюнайзена (табл. 8). Ошибка расчетов ограничивается точностью предположений о поведении тепловых свойств электронов.

Таблица 7.  

Энергия моделей никеля с потенциалом ЕАМ-B и с учетом электронных вкладов EeT

Т, К Z EeT, кДж/моль (6), (7)
1.0000 1.1111 1.1765 1.2500 1.3333 1.4286 1.5385 1.6667 1.8182
U + EeT
298 –420.57* –410.31* –396.56* –375.04* –343.52* –304.46* –239.94* –134.05* 36.01* 0.00
1000 –400.34* –390.34* –376.73* –355.35* –324.01* –284.72* –220.18* –114.36* 55.63* 2.06
2000 –366.57 –358.86* –345.97* –325.4* –294.73* –254.69* –190.29* –84.74* 84.94* 6.29
3000 –323.76 –311.47 –298.05 –276.94 –262.79* –220.51* –158.3* –53.17* 116.15* 12.10
5000 –261.31 –246.17 –231 –208.45 –175.95 –129.82 –63.52 43.1 186.03* 28.55
7000 –196.01 –178.71 –162.24 –138.45 –104.63 –56.88 10.42 118.18 289.44 51.93
10 000 –95.7 –75.57 –57.28 –31.89 3.05 52.55 121.82 231.21 404.74 93.23
15 000 79.4 103.24 123.92 151.86 189.44 240.57 344.71 424.74 600.83 176.14
20 000 266.13 294.48 315.73 345.63 384.91 438.41 513.15 627.43 805.95 275.31
25 000 468.73 498.45 521.86 553.51 595.14 650.35 726.82 843.73 1024.0 392.27

* Данные для устойчивых фаз ГЦК или ОЦК.

Таблица 8.  

Значения давления, модуля всестороннего сжатия и коэффициента Грюнайзена моделей никеля с потенциалом ЕАМ-B

Т, К Z
1.0000 1.1111 1.1765 1.2500 1.3333 1.4286 1.5385 1.6667 1.8182
p, ГПа
298 3.08* 32.35* 53.19* 79.82* 114.7* 160.7* 246.9* 406.3* 641.7*
1000 6.44* 34.88* 55.5* 82.26* 117.5* 164.7* 250.8* 409.8* 645.5*
2000 13.13 39.52* 59.85* 86.69* 122.6* 170.5* 256.5* 415.2* 651.6*
3000 23.86 50.57 70.56 97.38 128.3* 177.1* 262.8* 421.1* 658.1*
5000 33.76 61.52 82.22 109.57 146.3 196.8 284.6 442.5 672.5*
7000 41.86 70.66 91.95 119.9 157.3 209.2 298.4 456.5 695.2
10 000 52.47 82.71 104.9 133.8 171.9 225.8 317.2 475.7 715.8
15 000 67.75 99.93 123.5 153.9 193.9 250.4 344.7 504.5 746.6
20 000 81.21 115.8 140.1 171.8 213.3 272.7 369.7 530.6 774.9
25 000 93.92 129.7 155.3 188.3 231.7 293.7 392.9 555.0 801.2
KT, ГПа
298 200 377 422 456 564 872 1509 2433 2744
1000 188 371 420 461 575 881 1508 2421 2784
2000 177 263 416 494 593 869 1497 2455 2678
3000 192 311 391 498 608 864 1506 2425 2834
5000 181 343 410 489 628 914 1483 2522 2830
7000 194 366 428 495 631 939 1526 2396 2958
10 000 230 375 438 511 656 966 1545 2400 3005
15 000 228 392 460 540 692 1004 1573 2413 3055
20 000 263 404 473 563 728 1044 1600 2419 3132
25 000 249 418 493 592 763 1079 1623 2429 3184
γ
298 1.10 0.70 0.59 0.60 0.72 0.93 0.83 0.70 0.74
1000 1.10 0.81 0.72 0.71 0.79 0.92 0.83 0.71 0.74
2000 1.05 0.94 0.87 0.85 0.87 0.90 0.84 0.73 0.75
3000 1.01 0.89 0.87 0.86 0.94 0.88 0.84 0.75 0.75
5000 0.92 0.82 0.81 0.80 0.75 0.77 0.72 0.74 0.76
7000 0.84 0.76 0.75 0.74 0.71 0.73 0.69 0.70 0.67
10 000 0.73 0.68 0.67 0.67 0.64 0.67 0.65 0.64 0.62
15 000 0.55 0.55 0.55 0.55 0.55 0.57 0.58 0.56 0.54
20 000 0.40 0.43 0.43 0.44 0.46 0.49 0.52 0.48 0.47
25 000 0.26 0.33 0.33 0.34 0.38 0.41 0.45 0.40 0.40

* Данные для устойчивых фаз ГЦК или ОЦК.

В [66] проведено прямое моделирование процесса ударного сжатия атомной модели ГЦК-никеля при температурах до 1700 К с применением пяти различных потенциалов ЕАМ. Расчеты с потенциалом EAM-B дают для никеля при 1100 К и p = 0 значение γ = 1.46. Расчеты [66] с потенциалом [35] дают близкое значение γ = 1.37, а потенциал [26] и некоторые другие приводят к гораздо более высоким γ ≈ 2.0. Расчетные давления на адиабате в [66] существенно завышены против фактических при степени сжатия больше 1.3, а температуры примерно вдвое выше значений с потенциалами ЕАМ-A и ЕАМ-B.

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

Сравнение данных, полученных с потенциалами ЕАМ-А и ЕАМ-В, показывает, что роль тепловой энергии электронов весьма существенна. Например, при 298 К и Z = 1.8182 холодное давление ОЦК-никеля при учете электронных вкладов равно 641.7 ГПа, а без учета 583 ГПа (табл. 4). Температура на ударной адиабате значительно понижается при их учете (табл. 3). Наоборот, линия плавления мало чувствительна к наличию этих вкладов (табл. 6). Это позволяет приближенно рассчитывать линию плавления методом МД без учета электронных вкладов. Впервые проведены МД-расчеты ударных адиабат для твердых фаз и жидкого металла по отдельности, а также исследовано пересечение этих линий с линией плавления с учетом и без учета тепловых электронных вкладов. Плавление на ударной адиабате начинается при давлении ~276 ГПа (потенциал ЕАМ-В) или 272 ГПа (потенциал ЕАМ-А), т.е. здесь чувствительность к электронным вкладам также крайне низка. О правильности учета электронных вкладов можно судить по согласию с экспериментом расчетного холодного давления. В настоящем случае такое согласие (с [58]) есть (табл. 4).

В [58, 65] предполагается наличие в области высоких давлений только ГЦК-фазы никеля. Однако МД-анализ термодинамических свойств ОЦК, ГЦК и жидкого никеля показывает, что на фазовой диаграмме никеля при давлениях выше 100 ГПа вероятно появление области устойчивости ОЦК-фазы. Рассчитаны координаты тройной точки ГЦК–ОЦК–жидкость и предложена версия такой диаграммы (рис. 7).

Предложенный в данной работе потенциал ЕАМ должен достаточно точно описывать поведение никеля в интервалах температуры и давления, использованных при подборе параметров, т.е. при температурах до 25 000 К и степенях сжатия до 1.818. Для продвижения в область более высоких температур и давлений нужны дополнительные экспериментальные данные. Неясна (скорее, маловероятна) возможность применения этого потенциала для состояний с пониженной плотностью, когда вещество перестает быть металлом, в том числе в критической области и в виде газа.

Возможность использования данных ударного сжатия пористых образцов [73] подробно рассмотрена в [36], где показано, что сжатие пористого вещества не дает новой информации для состояний, покрываемых основной адиабатой, но расширяет интервалы достижимых параметров (T, p). Однако наличие пористости увеличивает вероятность получения неверной информации по чисто методическим причинам. Например, для никеля можно по известной основной адиабате практически точно рассчитать “теоретическую” ударную адиабату при пористости m = 1.108, но уже при m = 1.413 экспериментальная адиабата идет значительно выше теоретической [36]. Эти вопросы заслуживают специального анализа.

ЗАКЛЮЧЕНИЕ

Классический метод МД с применением потенциала ЕАМ позволил рассчитать термодинамические, структурные и диффузионные свойства твердого и жидкого никеля в широких интервалах температуры и давления, покрываемых основной ударной адиабатой. Учет тепловых вкладов электронов в энергию и давление сильно уменьшает температуру на адиабате и увеличивает холодное давление. Расчетная линия плавления моделей никеля плавно поднимается до 4518 К при 300 ГПа и почти нечувствительна к учету/неучету электронного вклада. Определены координаты участка плавления на ударной адиабате. Предложена фазовая диаграмма никеля с областью устойчивой ОЦК-фазы при давлениях выше 110–130 ГПа. Приведены таблицы энергии и давления моделей при степенях сжатия до Z = 1.8182, позволяющие рассчитать основные термодинамические свойства никеля. В отличие от непереходных металлов, точность расчетов ограничена из-за недостатка сведений о тепловых вкладах электронов в энергию и давление. Использование квантово-механических расчетов и данных ударного сжатия в МД-моделировании позволяет существенно продвинуться в исследовании свойств металлов в экстремальных условиях.

Список литературы

  1. Saito T., Sakuma Y. Densities of Pure Iron, Cobalt, and Nickel in the Molten State // J. Japan Inst. Metals. 1967. V. 31. № 10. P. 1140.

  2. Cheng Changrui, Xu Xianfan. Molecular Dynamics Calculation of Critical Point of Nickel // Int. J. Thermophys. 2007. V. 28. № 1. P. 9.

  3. www.webelements.com

  4. Gamsjäger H., Bugajski J., Gajda T. et al. Chemical Thermodynamics of Nickel. Nuclear Energy Agency Data Bank. V. 6. North Holland Publ., 2005.

  5. Desai P.D. Thermodynamic Properties of Nickel // Int. J. Thermophys. 1987. V. 8. № 6. P. 763.

  6. Mah A.D., Pankratz L.B. Contributions to the Data of Theoretical Metallurgy. XVI. Thermodynamic Properties of Nickel and its Inorganic Compounds. Washington: Bureau of Mines, 1976. Bull. 668.

  7. Chase M.W., Jr. NIST-JANAF Themochemical Tables. 4th ed. // J. Phys. Chem. Ref. Data. 1998. Monograph 9.

  8. Герасимов Я.И., Крестовников А.Н., Шахов А.С. Химическая термодинамика в цветной металлургии. Т. 4. М.: Металлургия, 1966. 427 с.

  9. Hixson R.S., Winkler M.A., Hodgdon M.L. Sound Speed and Thermophysical Properties of Liquid Iron and Nickel // Phys. Rev. B. 1990. V. 42. № 10. P. 6485.

  10. Abdullaev R.N., Kozlovskii Yu.M., Khairulin R.A., Stankus S.V. Density and Thermal Expansion of High Purity Nickel over the Temperature Range from 150 K to 2030 K // Int. J. Thermophysics. 2015. V. 36. № 4. P. 603.

  11. Assael M.J., Kalyva A.E., Antoniadis K.D. et al. Reference Data for the Density and Viscosity of Liquid Antimony, Bismuth, Lead, Nickel and Silver // High Temp.–High Press. 2012. V. 41. P. 161.

  12. Schroeder R.C., McMaster W.H. Compilation of Shock Hugoniot Data and Hugoniot Melting Points for the Elements // Lawrence Livermore Laboratory. UCRL-51 253.

  13. Филиппов С.И., Казаков Н.Б., Пронин Л.А. Скорость ультразвука, сжимаемость жидких металлов и их связь с различными физическими свойствами // Изв. вузов. Черная металлургия. 1966. № 3. С. 8.

  14. Waseda Y. The Structure of Non-Crystalline Materials: Liquids and Amorphous Solids. N.Y.: McGraw-Hill, 1980.

  15. http://res.tagen.tohoku.ac.jp/~waseda/scm/AXS/

  16. Eder O.J., Kunsch B., Suda M. et al. The Structure Factor of Liquid Nickel at 1873 K // J. Phys. F: Metal Physics. 1979. V. 9. № 7. P. 1215.

  17. Белащенко Д.К. Новый метод обращения структурного фактора некристаллической системы с помощью фурье-преобразования // Кристаллография. 1998. Т. 43. № 3. С. 400.

  18. LASL Shock Hugoniot Data / Ed. Marsh S.P. Berkeley: Univ. California Press, 1980.

  19. Al'tshuler L.V., Bakanova A.A., Dudoladov I.P. et al. Shock Adiabatic Curves of Metals // J. Appl. Mech. Tech. Phys. 1981. V. 22. P. 145.

  20. Al'tshuler L.V., Bakanova A.A., Trunin R.F. Shock Adiabats and Zero Isotherms of Seven Metals at High Pressures // Sov. Phys. JETP. 1962. V. 15. P. 65.

  21. http://www.ihed.ras.ru/rusbank

  22. Белащенко Д.К., Магидсон И.А. К расчету межчастичного потенциала в жидких никеле и железе в приближении Перкуса–Йевика // Изв. вузов. Черная металлургия. 1983. № 3. С. 4.

  23. Менделев М.И., Белащенко Д.К. Моделирование на ЭВМ структур жидкого Ni и аморфного сплава Ni62Nb38 // Металлы. 1995. № 3. С. 21.

  24. Landa A., Wynblatt P., Girshick A. et al. Development of Finnis–Sinclair Type Potentials for Pb, Pb–Bi, and Pb–Ni Systems: Application to Surface Segregation // Acta Mater. 1998. V. 46. № 9. P. 3027.

  25. Vitek V., Ackland G. J., Cserti J. // Materials Research Society. Pittsburgh. 1991. V. 186. P. 237.

  26. Mishin Y., Farkas D., Mehl M.J., Papaconstantopoulos D.A. Interatomic Potentials for Monoatomic Metals from Experimental Data and ab initio Calculations // Phys. Rev. B. 1999. V. 59. № 5. P. 3393.

  27. Baskes M.I. Determination of Modified Embedded Atom Method Parameters for Nickel // Mater. Chem. Phys. 1997. V. 50. P. 152.

  28. Zhang Q., Lai W.S., Liu B.X. Atomic Structure and Physical Properties of Ni–Nb Amorphous Alloys Determined by an n-body Potential // J. Non-Crystalline Solids. 2000. V. 261. P. 137.

  29. Koči L., Bringa E.M., Ivanov D.S. et al. Simulation of Shock-induced Melting of Ni Using Molecular Dynamics Coupled to a Two-temperature Model // Phys. Rev. B. 2006. V. 74. 012101.

  30. Özgen S., Sonğur L., Kara İ. Equations of State for Amorphous and Crystalline Nickel by Means of Molecular Dynamics Method // Turk J. Phys. 2012. V. 36. P. 59.

  31. Voter F., Chen S.P. Accurate Interatomic Potentials for Ni, Al, and Ni3Al // MRS Proc. 1986. V. 82. P. 175.

  32. Cherne F.J., Baskes M.I., Deymier P.A. Properties of Liquid Nickel: A Critical Comparison of EAM and MEAM Calculations // Phys. Rev. B. 2002. V. 65. 024209.

  33. Белащенко Д.К. Учет электронных вкладов в модели погруженного атома для железа и никеля при моделировании методом молекулярной динамики // ЖФХ. 2013. Т. 87. № 4. С. 633.

  34. Alemany M.M.G., Calleja M., Rey C. et al. A Theoretical and Computer Simulation Study of the Static Structure and Thermodynamic Properties of Liquid Transition Metals Using the Embedded Atom Model // J. Non-Cryst. Solids. 1999. V. 250–252. P. 53.

  35. Mendelev M.I., Kramer M.J., Hao S.G. et al. Development of Interatomic Potentials Appropriate for Simulation of Liquid and Glass Properties of NiZr2 Alloy // Phil. Magazine. 2012. V. 92. № 35. P. 4454.

  36. Belashchenko D.K. Liquid Metals. From Interatomic Potentials to the Properties, Shock Compression, Earth’s Core and Nanoparticles. Nova Science Publ., 2018.

  37. Daw M., Baskes M.I. Embedded-atom Method: Derivation and Application to Impurities, Surfaces, and other Defects in Metals // Phys. Rev. B. 1984. V. 29. № 12. P. 6443.

  38. Schommers W. Pair Potentials in Disordered Many-Particle Systems. A Study for liquid Gallium // Phys. Rev. A. 1983. V. 28. P. 3599.

  39. Doyama M., Kogure Y. Embedded Atom Potentials in fcc and bcc Metals // Comput. Mater. Sci. 1999. V. 14. P. 80.

  40. Girifalco L.A., Weizer V.G. Application of the Morse Potential to Cubic Metals // Phys. Rev. 1959. V. 114. P. 687.

  41. Ландау Л.Д., Лифшиц Е.М. Механика сплошных сред. М.: Гостехтеориздат, 1954. 796 с.

  42. Lin Zh., Zhigilei L.V. Electron-phonon Coupling and Electron Heat Capacity of Metals under Conditions of Strong Electron-phonon Nonequilibrium // Phys. Rev. B. 2008. V. 77. 075133.

  43. Lin Z., Zhigilei L.V. Temperature Dependences of the Electron-phonon Coupling, Electron Heat Capacity and Thermal Conductivity in Ni under Femtosecond Laser Irradiation // Appl. Surf. Sci. 2007. V. 253. P. 6295.

  44. Lord O.T., Wood I.G., Dobson D.P. et al. The Melting Curve of Ni to 1 Mbar // Earth Planetary Sci. Lett. 2014. V. 408. P. 226.

  45. Ландау Л.Д., Лифшиц Е.М. Статистическая физика. М.: ГИТТЛ, 1951.

  46. Al'tshuler L.V., Bushman A.V., Zhernokletov M.V. et al. // Sov. Phys. JETP. 1980. V. 51. № 2. P. 373.

  47. Levashov P.R., Sin’ko G.V., Smirnov N.A. et al. Pseudopotential and Full-electron DFT Calculations of Thermodynamic Properties of Electrons in Metals and Semiempirical Equations of State // J. Phys.: Condens. Matter. 2010. V. 22. P. 505501.

  48. Мигдал К.П. Термодинамические и кинетические свойства металлов с возбуждeнной электронной подсистемой. Автореф. дис. … канд. физ.-мат. наук. М., 2017.

  49. Simmons G., Wang H. Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook. Cambridge, Massachusetts: M.I.T. Press, 1971.

  50. Luo Fen, Chen Xiang-Rong, Cai Ling-Cang, Wu Qiang. Thermoelastic Properties of Nickel from Molecular Dynamic Simulations // J. At. Mol. Sci. 2011. V. 2. № 1. P. 10.

  51. Kamran S., Chen K.Y., Chen L. Ab initio Examination of Ductility Features of fcc Metals // Phys. Rev. B. 2009. V. 79. P. 024106.

  52. Chathoth S.M., Meyer A., Koza M.M., Juranyi F. Atomic Diffusion in Liquid Ni, NiP, PdNiP, and PdNiCuP Alloys // Appl. Phys. Lett. 2004. V. 85. P. 4881.

  53. Meyer A. The Measurement of Self-diffusion Coefficients in Liquid Metals with Quasielastic Neutron Scattering // EPJ Web of Conferences. 2015. V. 83. P. 01002.

  54. Jakse N., Wax J.F., Pasturel A. Transport Properties of Liquid Nickel near the Melting Point: An ab initio Molecular Dynamics Sudy // J. Chem. Phys. 2007. V. 126. 234 508.

  55. Cao Q.-L., Wang P.-P., Huang D.-H. et al. Properties of Liquid Nickel along Melting Lines under High Pressure // Chin. Phys. Lett. 2015. V. 32. № 8. P. 086201.

  56. Zharkov V.N., Kalinin V.A. Equations of State for Solids at High Pressures and Temperatures. N.Y.: Consultants Bureau, 1971.

  57. Urlin V.D. Melting at Ultra High Pressures in a Shockwave // Sov. Phys. JETP. 1966. V. 22. № 2. P. 341.

  58. Dewaele A., Torrent M., Loubeyre P., Mezouar M. Compression Curves of Transition Metals in the Mbar Range: Experiments and Projector Augmented-wave Calculations // Phys. Rev. B. 2008. V. 78. P. 104102.

  59. Dewaele A., Loubeyre P., Mezouar M. Equations of State of Six Metals above 94GPa // Phys. Rev. 2004. V. 70. 094112.

  60. Kormer S.B., Funtikov A.I., Urlin V.D., Kolesnikova A.N. Dynamic Compression of Porous Metals and the Equation of State with Variable Specific Heat at High Temperatures // Sov. Phys. JETP. 1962. V. 15. № 3. P. 477.

  61. McMahan K., Albers R.C. Insulating Nickel at 340 Million Atmospheres of Pressure // Phys. Rev. Lett. 1982. V. 49. P. 1198.

  62. Pozzo M., Alfè D. Melting Curve of Face-centered-cubic Nickel from First-principles Calculations // Phys. Rev. B. 2013. V. 88. № 2. P. 024111.

  63. Bouchet J., Bottin F., Jomard G., Zérah G. Melting Curve of Aluminum up to 300 GPa Obtained Through ab initio Molecular Dynamics Simulations // Phys. Rev. B. 2009. V. 80. 094102.

  64. Белащенко Д.К., Островский О.И. Молекулярно-динамическое моделирование ударного сжатия металлов. Железо и растворы железо‒сера // ЖФХ. 2011. Т. 85. № 6. С. 1063.

  65. Tonkov E.Yu., Ponyatovsky E.G. Phase Transformations of Elements under High Pressure. CRC Press, 2005. 377 p.

  66. Choi J., Yoo S., Song S. et al. Molecular Dynamics Study of Hugoniot Relation in Shocked Nickel Single Crystal // J. Mech. Sci. Technol. 2018. V. 32. № 7. P. 3273.

  67. Lazor P., Shen G., Saxena S.K. Laser-heated Diamond Anvil Cell Experiments at High Pressures: Melting Curve of Nickel up to 700 kbar // Phys. Chem. Minerals. 1993. V. 20. № 2. P. 86.

  68. Errandonea D. High-pressure Melting Curves of the Transition Metals Cu, Ni, Pd, and Pt // Phys. Rev. B. 2013. V. 87. P. 054108.

  69. Japel S., Schwager B., Boehler R., Ross M. Melting of Copper and Nickel at High Pressure: the Role of d-electrons // Phys. Rev. Lett. 2005. V. 95. P. 167801.

  70. Lord O.T., Wood I.G., Dobson D.P. et al. The Melting Curve of Ni to 1 Mbar // Earth Planetary Sci. Lett. 2014. V. 408. P. 226.

  71. Kerley G.I. Equations of State for Be, Ni, W, and Au // SANDIA Report. SAND 2003-3784. 2003.

  72. Minakov D.V., Levashov P.R. Melting Curves of Metals with Excited Electrons in the Quasiharmonic Approximation // Phys. Rev. B. 2015. V. 92. 224102.

  73. Трунин Р.Ф. Сжатие конденсированных веществ высокими давлениями ударных волн // УФН. 2001. Т. 171. С. 387.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Теплофизика высоких температур

Архивы выпусков Информация о журнале Отправить рукопись в журнал