ТЕПЛОФИЗИКА ВЫСОКИХ ТЕМПЕРАТУР, 2015, том 53, № 5, с. 683-691
УДК 536.4
ГИБРИДНЫЙ ПОТЕНЦИАЛ МЕЖЧАСТИЧНОГО ВЗАИМОДЕЙСТВИЯ И РАСЧЕТ ЛИНИЙ ПЛАВЛЕНИЯ ЛИТИЯ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ © 2015 г. Д. К. Белащенко
Национальный исследовательский технологический университет "Московский институт стали и сплавов" E-mail: dkbel@mail.ru Поступила в редакцию 21.02.2014 г.
Для металлических систем предложен многочастичный потенциал нового типа (ЕАМ-2), включающий в качестве параметра степень кристалличности атомов. Этот потенциал базируется на потенциале модели погруженного атома (ЕАМ). Степень кристалличности атома определяется путем разложения радиус-векторов соседних атомов по шаровым функциям и вычисления характеристики q6. Потенциал ЕАМ-2 хорошо описывает свойства жидкого, ОЦК- и ГЦК-лития при обычных и высоких давлениях. Рассчитаны линии равновесия ОЦК-литий—жидкость и ГЦК-литий—жидкость при давлениях до 40 ГПа и получено хорошее согласие с опытом, в том числе для максимума на кривой плавления. Приведены таблицы зависимости плотности и энергии моделей от давления и температуры. Обнаружено аномальное поведение плотности ОЦК- и ГЦК-лития в интервале давлений 25— 30 ГПа (рост плотности при изобарном нагреве).
Б01: 10.7868/80040364415040043
ВВЕДЕНИЕ
Расчет кривой равновесия кристалл—жидкость в широком интервале давлений методом молекулярной динамики — довольно непростая задача. Прямое определение температуры плавления в молекулярно-динамическом эксперименте затрудняется склонностью кристаллических моделей металла к перегреву. Расчет точки плавления путем вычисления энергий Гиббса жидкой и твердой фаз при различных температурах и нахождение точки пересечения этих графиков требует больших затрат времени. Однако основная трудность лежит в другом. При построении моделей жидкой и твердой фаз с одним и тем же межчастичным потенциалом соотношение плотностей этих фаз в точке плавления отличается от реального. Небольшая ошибка в равновесной плотности одной из фаз (например, всего на 1%) приводит к гораздо большей ошибке в изменении объема при плавлении АК Зависимость температуры плавления Тт от давления р описывается уравнением Клапейрона—Клаузиуса йТт/йр = ТАУ/АИ (АН — теплота плавления), так что ошибка в величине А V приводит к соответствующей ошибке наклона кривой равновесия йТт/йр, и даже если нормальная температура плавления рассчитана хорошо, то кривая равновесия при повышенных давлениях будет проведена неверно.
В случае металлических систем изменение объема при плавлении составляет обычно 1—3%
от всего объема, и получить здесь необходимую для расчета кривой плавления точность с одним и тем же межчастичным потенциалом практически невозможно. Уже отмечалось, что межчастичный потенциал, хорошо описывающий жидкую фазу, обычно недостаточно точен для твердой фазы [1]. Например, для жидкого лития свойства расплава вблизи точки плавления (Тт) довольно хорошо описываются потенциалом модели погруженного атома (ниже — потенциал ЕАМ-1 [2—5]). При 453.7 К и реальной плотности лития 0.5147 г/см3 [6] давление модели близко к нулю (0.00132 ± 0.0138 ГПа), а энергия составляет — 146.09 кДж/моль и близка к реальному значению —145.6 кДж/моль. Модуль всестороннего сжатия модели жидкого лития при нормальном давлении вблизи Тт равен 9.88 ГПа, что согласуется с данными по скорости звука (9.8 ГПа при 454 К [7]).
Однако ОЦК-литий при этих условиях описывается потенциалом ЕАМ-1 существенно хуже. Плотность реального ОЦК-лития при 298 К составляет 0.535 г/см3 [8], а плотность ОЦК-лития модели с потенциалом ЕАМ-1 при нулевом давлении равна 0.5514 г/см3, т.е. выше на 3%. Реальный модуль всестороннего сжатия лития равен 11 ГПа [8], а у модели при этих же условиях завышен (18.5 ГПа). При температуре 453.7 К с учетом коэффициента линейного расширения 46.1 х 10-6 К-1 получается плотность ОЦК-лития 0.5237 г/см3. Давление в модели ОЦК-лития при
Таблица 1. Расчет линии плавления лития. Потенциал ЕАМ-1 [5]
Давление, ГПа
0.00* 1.00 2.00 4.00 6.00 8.00 10.00
Жидкость
d 0.5147 0.5566 0.6084 0.6755 0.7295 0.7774 0.8213
U -146.09 -145.45 -145.07 -142.83 -139.66 -136.06 -132.14
ОЦК-литий
d 0.5382 0.5680 0.5973 0.6551 0.7099 0.7575 0.8010
U -150.56 -149.53 -148.92 -146.18 -142.64 -139.01 -134.98
Расчет по уравнению Клапейрона—Клаузиуса
AV 0.589 0.250 -0.212 -0.320 -0.263 -0.234 -0.214
AU 4.47 4.08 3.85 3.35 2.98 2.95 2.84
dTm/dp 59.8 30.6 -27.5 -47.8 -44.1 -39.7 -37.7
Примечания. d - плотность, г/см3; V - объем, см3/моль; p - давление, ГПа; A V - изменение объема при плавлении, см3/моль; AU - изменение энергии при плавлении, кДж/моль; dTm/dp - наклон линии плавления, К/ГПа. При расчетах производной dTm/dp принято ориентировочно T = 500 K. * При 453.7 K, AH = AU + pA V.
этих условиях равно —0.363 ГПа, т.е. модель заметно "растянута". При нулевом давлении равновесная плотность ОЦК-лития модели равна 0.5382 г/см3, т.е. на 2.7% выше фактической, а модуль всестороннего сжатия составляет 14.96 ГПа и также явно завышен. Энергия модели ОЦК-лития при 453.7 К равна—150.56 кДж/моль и теплота плавления АН = 150.56 - 146.09 = 4.47 кДж/моль, что заметно выше реального значения 3.0 кДж/моль [8]. Эти данные показаны в табл. 1.
Фактическое изменение объема при плавлении лития равно 6.94(0.5147-1—0.5237-1) = = 0.232 см3/моль. Изменение объема при плавлении модели ОЦК-лития с потенциалом работы [5] равно АУ = 0.589 см3/моль. По этим данным можно рассчитать производную йТт/йр по уравнению Клапейрона—Клаузиуса. Реальная величина производной равна 35—40 К/ГПа (см. [9]), а для модели лития йТт/йр = 59.8 К/ГПа. Поэтому увеличение давления на 1 ГПа приведет к росту расчетной температуры плавления на ~60 К вместо реального значения ~35 К. Следовательно, с ростом давления расхождения между расчетной линией двухфазного равновесия ОЦК—жидкость и фактической линией должны возрастать.
В табл. 1 приведены данные для расчета линии равновесия ОЦК-литий—жидкость, полученные путем моделирования методом молекулярной динамики (МД) [10] с потенциалом ЕАМ-1 работы [5]. Сжимаемость жидкости в моделях больше, чем у твердой фазы, и поэтому при давлениях выше 1 ГПа плотность жидкости выше плотности твердой фазы. По этой причине при давлениях 2—
10 ГПа производная йТт/йр всюду отрицательна, и с ростом давления температура плавления лития должна понижаться непрерывно и довольно быстро. Этот результат не согласуется с опытом, поскольку кривая плавления лития очень медленно поднимается с ростом давления [9].
Отсюда видно, что при описании свойств жидкой и кристаллической фаз одним и тем же потенциалом ЕАМ-1 не удается получить достаточную точность для описания линии равновесия кристалл—жидкость. Отмеченная здесь трудность хорошо известна и всегда возникает при обсуждении межчастичных потенциалов для ковалент-ных систем (углерод, кремний, германий и т.п.), где фазы с различной координацией атомов (связи sp3, зр2, ¿р и т.д.) описываются разными потенциалами. При исследовании таких систем, в которых координация атомов может быть различной, вводят многочастичные потенциалы, зависящие от локальной структуры, например от координационного числа [11, 12].
В настоящей работе впервые реализована аналогичная идея для металлов и сконструирован многочастичный потенциал ЕАМ-2 нового типа, который в случае жидкого лития должен совпасть с потенциалом ЕАМ-1, а в случае ОЦК-лития перейти в потенциал, хорошо описывающий твердую фазу. В качестве параметра, управляющего переходом из одной асимптотической формы потенциала в другую, выбрана "степень кристалличности" атомов системы.
Таблица 2. Величины 06 для различных структур [13, 14]
Структура Икосаэдр FCC НСР всс SC Жидкость
06 0.66332 0.57452 0.48476 0.51069 0.35355 0
МЕТОД РАСЧЕТА
Параметр порядка В качестве меры степени кристалличности использован параметр порядка 46, предложенный ранее в [13, 14]. Параметр порядка является комплексной величиной, определя-
емой как
1 Щ
Ят® = £^(0(г), Ф(г)).
у=1
(
(()
4п
21 +
1 £ \Ят ()|2
1
О, =
4п
21 +
п=-1
т=1
1 £ |°1т|
т=-1
ЯбтО
Ябт(®)
( т=6
£ |Ябт(')|2
V т=-6
у
Далее для пары (/ — соседний атом) определяется скалярное произведение векторов д6:
Здесь сумма берется по соседям атома /, находящимся в сфере заданного радиуса вокруг атома, а N - число этих соседей. В настоящих расчетах радиус сферы ближайших соседей для лития принят равным 4.00 А. Функции У1т(д, ф) — это шаровые функции, а полярный и азимутальный углы 9 и ф относятся к вектору г, соединяющему центральный и соседний атомы. При вычислении величин #т(0 в [13, 14] применяли сглаживание на границе сферы ближайших соседей. Усредняя локальные параметры д1т по всем частицам, получаем глобальный параметр 0,т:
N
О1т - ~ £ Я1т(®).
N
г=1
Далее можно сконструировать инварианты 47(0 и О,:
т=6
™ у = Яб(0 ' Яб(У) = £ Ябт(0Ябт * (У).
т=-6
Частицы / и]можно назвать "правильными соседями", если превышает некоторый порог (например, 0.5). Центральный атом является "твердоподобным", если число его правильных соседей превышает 6 [14]. Остальные атомы являются "жидкоподобными". Таким методом можно определять число жидкоподобных и твердопо-добных атомов в модели и наблюдать за процессами изменения структуры, кристаллизации или плавления. Этот метод был успешно применен, в частности, в работе [15] для исследования кластерной кристаллизации моделей серебра.
Потенциал модели погруженного атома. В настоящей работе для лития использовался потенциал модели погруженного атома ЕАМ-1, предложенный в [1—5]. В этом приближении потенциальная энергия системы записывается в виде
и = £ фр) + £ ф(Гу).
Величины 06 для различных структур показаны в табл. 2. Видно, что 06 мало зависит от структуры кристаллов, особенно при наличии тепловых колебаний решетки или дефектов, которые заметно понижают величину параметра порядка. Однако для жидкости этот параметр близок к нулю (~0.01). Поэтому 06вполне подходит для определения "степени кристалличности" системы [14]. Соответственно для определения степени кристалличности окружения отдельного а
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.