CC BY
64
УДК 538.95
Ю. С. Нагорнов, Р. Ю. Махмуд-Ахунов, М. Ю. Тихончев, Б. М. Костишко, В. Н. Голованов, В. В. Светухин
ПОСТРОЕНИЕ ТЕМПЕРАТУРНО-ЗАВИСИМОГО ПОТЕНЦИАЛА МЕЖЧАСТИЧНОГО ВЗАИМОДЕЙСТВИЯ ДЛЯ ДИОКСИДА УРАНА1
Аннотация. Применен новый подход к выбору потенциала межатомного взаимодействия для моделирования методом молекулярной динамики термических изменений свойств стехиометрического диоксида урана. В соответствии с теоремой Эренфеста и усредненным по времени потенциалом взаимодействия параметры потенциала в методе молекулярной динамики выбирались в виде медленно меняющихся температурных зависимостей. Рассчитанные значения постоянной решетки, энтальпии и теплоемкости при постоянном давлении с высокой точностью совпадают с экспериментальными данными в широком диапазоне температур от 250 до 3150 K, что подтверждает перспективность данного подхода.
Ключевые слова: молекулярная динамика, диоксид урана, потенциал взаимодействия.
Abstract. The new approach was applied to the construction of the atom interaction potential for molecular dynamics modeling of thermal changes of properties of uranium dioxide. The potential parameters were selected in the form of slowly varying temperature dependences in accordance with Ehrenfest theorems and average potential. This method was able to calculate of the temperature dependences of constant lattice, enthalpy and heat capacities up to 3150K with hi accuracy. These results are confirmed perspective of presented approach.
Keywords: molecular dynamics, uranium dioxide, interaction potential.
Введение
Метод молекулярной динамики (МД) широко используется для моделирования поведения материалов под облучением [1], поскольку получение экспериментальных данных сопряжено с большими затратами. Стехиометрический диоксид урана является наиболее используемым видом ядерного топлива в промышленных реакторах и достаточно хорошо экспериментально изученным актинид-оксидным топливом [2]. В настоящее время активно ведутся исследования свойств этого материала с помощью метода МД. Так, в работе [3] показано, что процессы миграции и кластеризации вакансий атомов кислорода протекают наиболее активно вблизи границы дислокаций, что является одним из возможных путей формирования микроструктуры оболочки таблетки топлива при выгорании в реакторе. Кроме этого, метод МД позволяет исследовать радиационно-стимулированную диффузию [4], динамику фазовых переходов под давлением [5] и многое другое.
Существенным ограничением МД моделирования является ограничение моделируемого времени поведения системы, не превышающего в лучшем
1 Работа выполнена в рамках реализации ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг. и при поддержке гранта РНПВШ 2.1.2/5656.
случае сотен наносекунд. Однако этот недостаток можно преодолеть при исследовании миграционных процессов с помощью приема, называемого методом термически ускоренной динамики [3]. Зная энергию активации процесса, можно перенести данные расчета, полученные при высоких температурах, на низкие температуры с одновременным увеличением моделируемого времени на несколько порядков. В связи с этим расчеты при высоких температурах становятся особенно актуальными, поскольку позволяют проводить моделирование методом МД на времена до десятков микросекунд.
Для моделирования свойств диоксида урана методом МД используется парный потенциал межатомного взаимодействия с параметрами, определяемыми, как правило, полуэмпирическим методом [6]. Выбор потенциала взаимодействия частиц является наиболее важным этапом построения численной модели и, как правило, вызывает немало дискуссий. Вообще говоря, для нахождения потенциала взаимодействия необходимо решить уравнение Шре-дингера для рассматриваемой системы частиц. Однако, как показала практика, использование потенциалов, получаемых из первых принципов, для МД расчетов приводит к отрицательным результатам. Поэтому широкое распространение получил другой подход, связанный с использованием простых модельных потенциалов, вид которых, как правило, основывается на теоретических представлениях о наиболее существенных вкладах в рассматриваемый тип взаимодействия. Параметры потенциалов находятся из сопоставления расчетных и экспериментальных данных. При этом потенциал, оптимизированный по одному свойству, может существенно отличаться от потенциала, оптимизированного по другому свойству. То есть выбор потенциала обосновывается конкретной задачей.
С развитием вычислительной техники потенциалы становятся более сложными, а количество параметров постоянно увеличивается. Например, в работе [7] при моделировании магнитных свойств в потенциале появляются слагаемые, учитывающие электронную плотность для различных спиновых состояний, что должно привести к температурной зависимости параметров потенциала. В общем случае потенциал взаимодействия можно определить, зная потенциальную энергию атома в каждой точке кристалла. В соответствии с теоремой Эренфеста для среднего состояния механических величин в одномерном случае имеет место квантовое уравнение Ньютона [8]:
d 2 х дU (x)
М----=-------—,
dt дx
где М - центр масс волнового пакета атома; U (x) - усредненная по времени
потенциальная энергия атома в кристалле; x - координата центра тяжести волнового пакета атома в кристалле.
В общем случае после применения правила вычисления среднего значения силы в квантовой механике и разложения в ряд Тейлора по степеням Лx можно получить [8]
дЦЫ =ЭU(X) +1 дЬ(^ +
дx дx 2 дк3
При моделировании методом МД принято пользоваться классическими уравнениями Ньютона, где учитывается только первое слагаемое усреднен-
й 1 d3U (X) —2
ной силы, в то время как малая величина-----------—— также дает вклад
2 Эх3
в динамику движения. Необходимо заметить, что первое слагаемое будет неизменным, в то время как второе будет расти с ростом температуры и дисперсии колебаний атомов, и его вклад будет увеличиваться. Именно поэтому в настоящей работе для моделирования термических изменений свойств стехиометрического диоксида урана предлагается использовать потенциал с параметрами, зависящими от температуры. Аналогичный подход был использован авторами работ [8, 9], в которых применялся эффективный потенциал взаимодействия при моделировании сплава AuCd.
1. Метод моделирования
Моделирование методом МД проводилось в периодических граничных условиях с использованием программного комплекса Moldy [10]. Транслируемая ячейка была выбрана в виде кубического кристалла со структурой флюорита и содержала 768 ионов, или 4*4*4 элементарных ячеек. Для всех расчетов шаг интегрирования был выбран 4 фс, радиус обрезания потенциала был равен 1 nm. В зависимости от задачи расчеты велись либо в ансамбле NVE при постоянном объеме и энергии, либо в ансамбле NPT при постоянном давлении и температуре, количество частиц в системе было неизменной величиной в обоих случаях.
Потенциал межатомного взаимодействия был выбран в форме потенциала Борна - Майера, что обеспечило минимальный набор параметров, некоторые из которых взяты в виде кусочно-линейных функций температуры:
2 АЛ
ai + aj ~ rij
bi + bj
Zi (T)Z: (T)e
U (rij) = ------------j---+ f0(T )(bi + bj )exp
rj
cicj
r6 ij
где первое слагаемое соответствует кулоновскому взаимодействию, а второе и третье - потенциалу Борна - Майера [12]. Значения не зависимых от температуры параметров потенциала были взяты из работы [13].
Для восстановления параметров потенциала Zl и /0 были использованы экспериментальные данные по тепловому расширению решетки Ь02 и изменению энтальпии [11]. Интересно, что полученные по этой методике значения параметров дробного заряда и/0 имеют две характерные линейные области с переходом при температуре вблизи 2500 К температуре перехода в су-перионное состояние (рис. 1) [14]. Подобранные значения параметров потенциала аппроксимировались линейными зависимостями в каждой области (до и после 2500 К). Во всех дальнейших расчетах значения параметров вычислялись из полученных кусочно-линейных зависимостей. Подобная методика введения линейных температурных зависимостей позволяет существенно упростить процедуру подбора параметров потенциала для других материалов. В подтверждение предлагаемой методики расчета можно заметить, что температурное изменение параметров потенциала существенно ниже их абсолютных величин, т.е. проведенный расчет показывает правильность принятого приближения при выборе потенциала.
0,З8 0,З7 -І
Е 0.З6Н
0,Зб-0,З4 -0,ЗЗ
1500
2000
2б00
З000
ф
03
о
о
ф
О)
о
T, K
Рис. 1. Зависимость параметраf и дробного заряда иона кислорода от температуры и кусочно-линейная аппроксимация с переходом вблизи точки 2500 K
2. Расчетные оценки энтальпии и термического расширения кристалла
Предложенная в настоящей работе методика выбора потенциала (далее потенциал NMKG) позволяет рассчитать параметр решетки и энтальпию диоксида урана с достаточной точностью во всем температурном интервале (рис. 2, 3). При этом параметры потенциала NMKG меняются не более чем на 10 % (рис. 1). Необходимо отметить, что зависимости zt и f0, рассчитанные из наиболее применяемых на сегодня потенциалов, таких как Basak [15], Morelon [16], Yamada [17], существенно отличаются от экспериментальных в диапазоне температур 1500-3150 K, тем более очевидны отличия для классических потенциалов Arima и Lewis [11]. Поскольку основные расчеты велись для температурного интервала 600-1100 K, этими отклонениями до сих пор пренебрегали. Однако с использованием метода термически ускоренной динамики температурный интервал может быть существенно расширен в область высоких температур.
Интересно, что параметр решетки, рассчитанный с использованием потенциалов Basak и Morelon, совпадает с экспериментальными данными достаточно хорошо (рис. 2), но при этом функция энтальпии расходится с ростом температуры (рис. 3). Напротив, значения, полученные с использованием потенциала Yamada, дают существенную погрешность при расчете параметра решетки, но при подборе модуля упругости приводят к наиболее точному результату.
Отличие численных и экспериментальных значений вынуждает искать новые подходы, что часто приводит к необоснованному увеличению количества параметров в потенциале [1, 11]. Сегодня только форм потенциалов известно более шести [10], что, соответственно, приводит к различиям полученных на их основе численных результатов и свидетельствует о неоднозначности метода МД. Предлагаемый в настоящей работе подход позволяет
і59
надеяться, что учет температурных зависимостей параметров потенциала, независимо от его выбранной формы, должен приводить к одинаковым численным результатам и точному соответствию известным экспериментальным данным.
5,70 < 5,65-
Ф
-55 5,60
Е га
го 5,55
ш
о
го
5,50
5,45
5,40
-- а NMKG Experime х Basak + Morelon □ Yamada д Arima ■ Lewis
А
А
А'
>
□
\/| У □ ПА А
и А А
а й Л
і ^ А ■
'Л1 ■ ■
■ ■
■ ■
■ ■
і
500 1000 1500 2000 2500 3000 3500
Т, К
Рис. 2. Сравнение экспериментальных данных температурной зависимости параметра решетки и02 с данными расчета настоящей работы (КМКв) с использованием различных потенциалов [15-17]
300
о 250 Е
£ 200 го
с
Ш
150
100
■ N№1^ Experiment Д Basak а Potashnikov о МогеІоп ■
/ /
/ А А
✓ У л * с 6 О
✓ у г 6^°
✓ ✓ Ж / •> о
\ V \ ш \ ''ь°с
.Л'\ А
2 1 1
1400 1600 1800 2000 2200 2400 2600 2800 3000 3200
Т, К
Рис. 3. Сравнение экспериментальных данных температурной зависимости энтальпии с данными расчета настоящей работы (КМКв) с использованием различных потенциалов [15, 16, 18]
Предложенный в настоящей работе подход дает хорошее согласие с экспериментальными значениями во всем диапазоне температур, при этом погрешность расчетов не превышает 0,5 % (рис. 3). Погрешность расчетов с применением наиболее точных на сегодня потенциалов варьируется от 2 до 17 % в зависимости от температуры в диапазоне 1500-3000 К. Необходимо отметить, что получено хорошее согласие двух экспериментальных зависимостей с расчетами на основе потенциала КМКО. Для проверки подбора параметров /0 и г в настоящей работе были проведены дополнительные расчеты теплоемкостей СР и Су.
3. Расчет теплоемкостей Ср и Су
Теплоемкости при постоянном давлении и различных температурах рассчитывали, моделируя ансамбль КРТ. Полученные значения энергии системы атомов аппроксимировались функцией от температуры в виде полинома. Затем производная от энергии по температуре вычислялась как производная от полученных полиномов. Как показали расчеты (рис. 4), зависимости, вычисленные с использованием потенциала КМКО, также дают хорошее согласие с экспериментальными значениями во всем диапазоне температур, в отличие от потенциалов других авторов [11], где погрешность расчетов варьируется в широком интервале 10-90 %.
180 -160140-
K
*ol 12G Е
“ -100-BG 6G 4G
О
1 1
о Experiment х Basak ■ Morelon + Yamada a NMKG о Arima Л Lewis ▲
і
k
/
і Y
nfS F
■
+ ■ V* н*Л. * ¥ *
XO&X С Lfee* <■ <^XB ■ X* <r с с Л о ло z о 4с
°аа гА
7
і
5GG
1000
1500 T, К
2GGG
25GG
3000
Рис. 4. Зависимость теплоемкости от температуры, вычисленная при постоянном давлении в ансамбле ЫРТ с использованием различных потенциалов [15-17]. КМКв - данные расчета настоящей работы
Необходимо отметить, что в расчетах с более простыми потенциалами, таким как Arima и Lewis [11], получается теплоемкость, линейно зависимая от температуры и меняющаяся в диапазоне 75-85 J/(mol*K). Искусственное добавление полиномов 3-5 степени в форму потенциала не отвечает какому-
l6l
либо физическому процессу, но позволяет получать небольшой перегиб в области 2000 К [11, 16, 17]. Предложенный подход с потенциалом NMKG позволил получить не только соответствующую экспериментальной форму температурных зависимостей теплоемкостей CP и CV, но и рассчитать их численные значения с высокой точностью. При этом, несмотря на то, что параметры потенциала подбирались в ансамбле NVE, более высокая точность расчетов была получена в ансамбле NPT для теплоемкости CP.
Заключение
В работе предложен подход температурной зависимости параметров потенциала межатомного взаимодействия при моделировании свойств диоксида урана. Показано, что этот подход позволяет существенно повысить точность расчетов в методе молекулярной динамики. С использованием предлагаемого подхода и потенциала Борна - Майера были рассчитаны температурные зависимости постоянной решетки, энтальпии и теплоемкости при постоянном давлении. Получено хорошее согласие с экспериментальными данными, в том числе при расчетах в высокотемпературной области, где результаты других исследователей существенно отклоняются от эксперимента.
Авторы благодарят профессора Ульяновского государственного университета В. М. Журавлева за ценные замечания и обсуждение результатов работы.
Список литературы
1. Tikhonchev, M. MD simulation of atomic displacement cascades in Fe-10 at.%Cr binary alloy / M. Tikhonchev, V. Svetukhin, A. Kadochkin, E. Gaganidze // J. Nuclear Materials. - 2009. - V 395. - P. 50-57.
2. Маершин, А. Тепловыделяющие элементы с виброуплотненным оксидным топливом / А. Маершин. - Димитровград : ФГУП ГНЦ РФ НИИАР, 2007. - 327 с.
3. Ichinomiya, T. Temperature accelerated dynamics study of migration process of oxygen defects in UO2 / T. Ichinomiya, B. Uberuaga, K. Sickafus // J. Nuclear Materials. - 2009. - V. 384. - P. 315-321.
4. Martin, G. A molecular dynamics study of radiation induced diffusion in uranium dioxide / G. Martin, S. Maillard, L. Brutzel, P. Garcia, B. Dorado, C. Valot // J. Nuclear Materials. - 2009. - V. 385. - P. 351-357.
5. Desai, T. G. Molecular dynamics study of diffusional creep in nanocrystalline UO2 / T. G. Desai, B. P. Uberuaga, D. Wolf // Acta Materialia. - 2008. - № 56. - P. 44894497.
6. Методы молекулярной динамики в физической химии / под ред. Ю. М. Товбина. -М. : Наука, 1996. - 334 с.
7. Ackland, G. Two-band second moment model for transition metals and alloys / G. Ackland // J. Nuclear Materials. - 2006. - V. 351. - Р. 20-27.
8. Блохинцев, Д. И. Основы квантовой механики / Д. И. Блохинцев. - М. : Наука, 1976. - 664 с.
9. Guthikonda, V. S. An Effective Interaction Potential Model for the Shape Memory Alloy AuCd / V. S. Guthikonda, R. S. Elliott // Continuum Mechanics and Thermodynamics. - 2009. - V. 21(4). - Р. 269-295.
10. Govers, K. Comparison of interatomic potentials for UO2. Part I: Static calculations / K. Govers, S. Lemehov, M. Hou, M. Verwerft // Journal of Nuclear Materials. - 2007. -V. 366. - P. 161-177.
11. Govers, K. Comparison of interatomic potentials for UO2 Part II : Molecular dynamics simulations / K. Govers, S. Lemehov, M. Hou, M. Verwerft // Journal of Nuclear Materials. - 2008. - V. 376. - P. 66-77.
12. Гулд, Х. Компьютерное моделирование в физике / Х. Гулд, Я. Табочник. - М. : Мир, 1990. - Ч. 1. - 349 с.
13. Yamasaki, S. Evalution of Thermal Conductivity Hyperstoihiometric UO2+x by Molecular Dynamics Simulation / S. Yamasaki, T. Arima, K. Idemitsu [et al.] // International Journal of Thermophysics. - 2007. - V. 28. - № 2. - P. 661-673.
14. Молодец, А. М. Фазовые переходы диоксида урана при высоких температурах и давлении / А. М. Молодец, В. Е. Фортов // Письма в ЖЭТФ. - 2004. - Т. 80. -№ 3. - С. 196-199.
15. Basak, C. B. Classical molecular dynamics simulation of UO2 to predict thermophysical properties / C. B. Basak, A. K. Sengupta, H. S. Kamath // J. Alloys and Comp. - 2003. - V. 360. - P. 210-216.
16. Morelon, N.-D. A new empirical potential for simulating the formation of defects and their mobility in uranium dioxide / N.-D. Morelon, D. Ghaleb // Phil. Mag. - 2003. -V. 83. - P. 1533-1550.
17. Yamada, K. Evaluation of thermal properties of uranium dioxide by molecular dynamics / K. Yamada, K. Kurosaki, M. Uno, S. Yamanaka // J. Alloys and Comp. -2000. - V. 307. - P. 10-15.
18. Поташников, С. И. Молекулярно-динамическое восстановление межчастич-ных потенциалов в диоксиде урана по тепловому расширению / С. И. Поташников, А. С. Боярченков, К. А. Некрасов, А. Я. Купряжкин // Альтернативная энергетика и экология. - 2007. - № 8 (52). - C. 43-52.
Нагорнов Юрий Сергеевич кандидат физико-математических наук, доцент, начальник лаборатории, Ульяновский государственный университет
Е-таі1: nagomovys@ulsu.ru
Махмуд-Ахунов Руслан Юсупович младший науный сотрудник, Ульяновский государственный университет
Е-таі1: nagornovys@u1su.ru
Тихончев Михаил Юрьевич
кандидат физико-математических наук, начальник лаборатории, Ульяновский государственный университет
Е-таі1: tikhonchev@sv.u1su.ru
Костишко Борис Михайлович
доктор физико-математических наук, профессор, ректор, Ульяновский государственный университет
Е-таі1: contact@u1su.ru
Nagornov Yuriy Sergeevich Candidate of physical and mathematical sciences, associate professor, head of laboratory, Ulyanovsk State University
Makhmud-Akhunov Ruslan Yusupovich Junior researcher, Ulyanovsk State University
Tikhonchev Mikhail Yurievich Candidate of physical and mathematical sciences, head of laboratory,
Ulyanovsk State University
Kostishko Boris Mikhaylovich Doctor of physical and mathematical sciences, professor, rector of Ulyanovsk State University
1б3
Ульяновский государственный университет
Голованов Виктор Николаевич
доктор физико-математических наук, профессор, проректор по науке и информационным технологиям,
Golovanov Viktor Nikolaevich
Doctor of physical and mathematical sciences, professor, vice rector for scientific affairs and information technologies, Ulyanovsk State University
E-mail: golovanovvn@ulsu.ru
Светухин Вячеслав Викторович
доктор физико-математических наук, профессор, директор научноисследовательского технологического института Ульяновского
Svetukhin Vyacheslav Viktorovich
Doctor of physical and mathematical sciences, professor, director of the research technological institute under Ulyanovsk State University
государственного университета E-mail: slava@sv.uven.ru
УДК 538.95 Нагорнов, Ю. С.
Построение температурно-зависимого потенциала межчастичного взаимодействия для диоксида урана / Ю. С. Нагорнов, Р. Ю. Махмуд-Ахунов, М. Ю. Тихончев, Б. М. Костишко, В. Н. Голованов, В. В. Светухин // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2010. - № 3 (15). - С. 156-164.
