CC BY
62
УДК 621.039531
МОДЕЛИРОВАНИЕ ГПУ-ЦИРКОНИЯ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
© 2013 П.Е. Капустин Ульяновский государственный университет Поступила в редакцию 26.11.2013
Работа направлена на изучение радиационного повреждения ГПУ-7г методом молекулярной динамики. Проводится расчет энергии формирования точечных дефектов, а также пороговой энергии смещения. Моделируется каскад атомных смещений при различных температурах, энергиях и потенциалах межатомного взаимодействия. Полученные результаты сопоставляются с имеющимися литературными данными.
Ключевые слова: цирконий, метод молекулярной динамики, точечные дефекты, пороговая энергия смещения, каскад атомных смещений.
ВВЕДЕНИЕ
Атомная энергетика активно развивается в последние десятилетия. Совершенствование технологий, конструкции и безопасности ядерных реакторов необходимо для дальнейшего ее развития. Однако без глубокого и детального изучения взаимодействия излучения с веществом это невозможно. Данные исследования позволят не только расширить представления о физической природе процессов, происходящих в материалах под облучением, но и предоставит возможности использования новых материалов и методов при конструировании ядерных реакторов, а также модернизировать уже имеющиеся реакторные установки.
Развитие вычислительных средств в последние годы позволяет рассматривать компьютерное моделирование как наиболее удобный и перспективный метод исследований в области реакторного материаловедения [1]. Вычислительный эксперимент позволяет детально изучать различные процессы, происходящие в материалах под облучением, в то время как проведение натурных экспериментов зачастую затруднено, а в ряде случаев - невозможно. Компьютерное моделирование основывается на теоретических основах происходящих процессов и позволяет изучать не только равновесное положение, но и динамику протекающих таких процессов.
Одним из основных методов вычислительного эксперимента, позволяющих исследовать процессы на уровне отдельных атомов, является метод молекулярной динамики (МД). Этот метод основывается на решении уравнений механики Ньютона для совокупности взаимодействующих атомов, образующих модельный кристаллит.
Капустин Павел Евгеньевич, аспирант, младший научный сотрудник Научно-исследовательского технологического института. E-mail: [email protected]
Настоящая работа посвящена МД-моделиро-ванию ГПУ циркония. Цирконий обладает удачным сочетанием ядерно-физических характеристик и механических свойств и потому считается одним из наиболее перспективных конструкционных материалов ядерных энергетических установок. В настоящее время сплавы циркония широко используются для изготовления оболочек топливных стержней и элементов конструкции тепловыделяющих сборок реакторов с водой под давлением.
В работе приводятся результаты расчетов энергий формирования дефектов кристаллической структуры, оценивается пороговая энергия смещения, моделируются каскады атомных смещений и подсчитываются количества пар Френкеля, переживших рекомбинацию в каскаде («выжившие» дефекты). Полученные в ходе расчетов результаты сопоставляются с имеющимися в литературе данными.
МЕТОД МОДЕЛИРОВАНИЯ
Для циркония в данной работе используется многотельный потенциал, предложенный Мен-делевым и Акландом в работах [2, 3]. В работе [3] авторами предлагается несколько потенциалов. Обозначенный ими как потенциал #3 используется для ГПУ-2г, который и рассматривается в данной работе. В литературе отмечается [4, 5, 6], что данный потенциал является наиболее совершенным и дающим хорошее согласие с экспериментальными данными. Использование потенциала, приведенного в работе [2] сделано с целью сопоставления различия результатов моделирования и возможного выявления особенностей или недостатков модернизированного потенциала по отношению к его более ранней версии. Оба являются потенциалами вида Финни-са-Синклера [7].
Для проведения дальнейших расчетов использовалась модель кристаллита ГПУ-2г, а также периодические граничные условия. Построение модели осуществлялось с помощью трансляционной симметрии, задав изначально координаты базовых точек. В зависимости от температуры варьировались параметры решетки.
ЭНЕРГИЯ ФОРМИРОВАНИЯ ТОЧЕЧНЫХ ДЕФЕКТОВ
При воздействии различных факторов (например, нагрев, облучение и т.п.), в кристаллитах могут образовываться точечные дефекты (нарушение трансляционной симметрии кристалла): вакансии и собственные междоузельные атомы (СМА). Используя методы компьютерного моделирования можно добавить или убрать один из атомов. Реализовать такое на практике не представляется возможным, однако для исследовательских целей такой метод не только прост, но и достаточно информативен.
В таблице 1 приведены результаты вычислительных экспериментов, описанных в работах [2, 8, 9]. Рассмотренные конфигурации СМА в ГПУ решетке циркония показаны на рис. 1.
Отсутствие данных по той или иной конфигурации объясняется тем, что такая конфигурация неустойчива и при релаксации она переходит в
другую рассматриваемую конфигурацию, что так же отмечалось в указанных выше работах
ПОРОГОВАЯ ЭНЕРГИЯ СМЕЩЕНИЯ
При взаимодействии с материалом происходят упругие и неупругие столкновения бомбардирующих частиц с атомами. В результате таких столкновений часть энергии излучения передается атомам, которые называются первично выбитыми атомами (ПВА). Если переданная энергия оказывается достаточно большой, атом может сместиться из своего узла кристаллической решетки и образовать пару Френкеля (вакансия - междоузельный атом). Минимальная кинетическая энергия, необходимая для образования устойчивой пары Френкеля, называется пороговой энергией смещения.
В кристаллите выбирается ПВА, которому сообщается определенная энергия и направление движения. Устойчивой парой Френкеля в рамках данной работы считалась такая пара, которая не прорекомбинировала в течение 10 пикосекунд. Расчет проводился для четырех изотропных случайных направлений на двух потенциалах [2, 3]. Как было сказано выше, из работы [3] брался потенциал для ГПУ-2г (#3). Для каждого рассматриваемого направления ПВА строился кристаллит своих размеров для избегания «зацикливания» возмущений вследствие периодических граничных условий.
Конфигурация [2] [8] [9]
Вакансия 1.786 - -
ВС 3.756 3.24 -
В8 3.760 3.53 3.60
ВО 3.970 3.24 3.53
С 3.979 3.56 -
8 4.319 3.44 4.57
Б С
ВО ВЭ ВС
Рис. 1. Конфигурации собственного СМА в ГПУ решетке циркония
Таблица 1. Энергии формирования вакансии и различных конфигураций СМА в ГПУ решетке циркония, эВ
Количество атомов в кристаллите варьировалось от 12 до 24 тысяч в зависимости от условий вычислительного эксперимента.
Для нахождения пороговой энергии проводилась серия экспериментов с шаговым увеличением энергии ПВА на 1 эВ. Как только удавалось получить пару Френкеля, которая не прореком-бинировала в течение 10 пикосекунд, данная энергия ПВА фиксировалась. Ввиду того, что даже превысив пороговую энергию, можно получить ситуацию, где пара Френекеля прорекомбениро-вала быстрее, шаговый метод является более достоверным, исключая ошибку перескока. Резуль-
таты моделирования приведены в табл. 2.
Пороговая энергия сильно зависит от выбранного кристаллографического направления, что определяется кристаллическим строением.
КАСКАДЫ СМЕЩЕНИЙ
При сообщении некоторому атому достаточно большой энергии инициируется каскад атомных смещений. Размер каскадной области определяется энергией ПВА. На рис. 2 представлен пример эволюции дефектной структуры кристаллита 2г при похождении каскада с энергией 5 кэВ.
Таблица 2. Пороговая энергия смещения
Направление ПВА <0001> < 1210 > <1100> <0П0>
Потенциал [2] [3] [2] [3] [2] [3] [2] [3]
Пороговая Е, еУ 50 63 34 34 65 77 119 124
Рис. 2. Эволюция дефектной структуры кристаллита при развитии каскада смещений: а - 0,09 рв; Ь - 0.26 пс; с - 0,65 пс; а - 3,85 пс; е - 7,52 пс; £ - 44.84 пс Синие - вакансии, черные - междоузельные атомы
Пик каскада, завещающий так называемую баллистическую стадию, приходится на 0,65 пи-косекунд. Миновав его, начинается процесс ре-комбенации точечных дефектов, по окончанию которого остается относительно небольшое количество дефектов («выжившие»дефекты).
В ходе экспериментов в данной работе использовались периодические граничные условия. Исходя из этого размеры кристаллита, а также энергия ПВА подбиралась так, чтобы не возник эффект наложения каскада на самого себя, что означало бы получение заведомо неверных данных.
Случайным образом определялось 8 направлений для ПВА. Расчет проводился для 1, 5, 10, 15 и 20 кэВ с размерами кристаллита ~100, 200, 300, 400, 500 тысяч атомов соответственно. Для расчета использовался потенциал из работы [3]. Результаты представлены в таблице 3, где р и О это углы в сферической системе координат, определяющие направление движения ПВА: р -угол, отсчитываемый от направления <1210 > в плоскости (0001), О - угол, отсчитываемый от направления <0001>, Е - энергия ПВА, N - количество атомов в модели кристаллита, Ь -моделируемое время каскада. Перед моделированием каскада кристалл релаксировали при тем-
пературе 600 К и нулевом давлении в течение 10 пс. Указанные в табл. 3 и далее погрешности являются статистическими погрешностями, соответствующими одному стандартному отклонению для выборочного среднего (1 о ). Также был произведен расчет эффективности каскада по NRT-стандартy
Также были проведены расчеты с помощью потенциала из работы [2] для энергий ПВА 15 и 20 кэВ соответственно при той же температуре. Полученные результаты представлены в табл. 4.
Полученные результаты показывают, что из общего числа атомных смещений от ПВА с энергией от 10 до 20 кэВ выживает 15-18% дефектов при работе с потенциалом [3] и 18-20% при использовании потенциала [2].
Изучения влияния температуры на количество выживших дефектов проводилось в работах [4, 10]. Отмечалось, что с ее ростом количество дефектов не изменялось в пределах погрешности. Проведенные исследования в рамках данной работы в диапазоне температур от 100 до 400 К согласуются с полученными ранее данными. Соответствующие результаты приведены в табл. 5. Для расчетов использовался потенциал из работы [3], и энергия ПВА 10 кеУ. В табл. 6 указаны
Таблица 3. Количество выживших дефектов
Направление Е = 1 кеУ № 104832 г=23 РБ Е = 5 кеУ N=209880 1=36 РБ Е= 10кеУ N=296888 1=45 РБ Е=15 кеУ N=405080 1=70 РБ Е=20кеУ N=492800 1=88 рб
ф = 43.314 0 = 70.999 4 9 4 27 29
ф = 74.023 0 = 71.563 4 13 18 40 25
ф = 48.309 0 = 59.019 4 6 16 27 54
ф = 61.289 0 = 30.179 3 6 10 9 40
ф = 67.271 0 = 82.287 1 2 30 22 71
ф = 60.42 0 = 54.605 3 6 26 47 3
ф = 66.783 0 = 62.429 3 6 11 24 9
ф = 50.196 0 = 41.611 2 7 8 21 35
Среднее 3±0.378 6.875±1,109 15,375±3,179 27,125±4,147 33,25±7.898
Р 0.3±0.038 0.138±0.022 0.154±0.032 0.181±0.028 0.166±0.039
Таблица 4. Количество выживших дефектов
Направление E = 15 keV, N = 405080, t=70 ps E = 20 keV, N = 492800, t = 88 ps
ф = 43.314 0 = 70.999 39 41
ф = 74.023 0 = 71.563 29 34
ф = 48.309 0 = 59.019 25 35
ф = 61.289 0 = 30.179 35 42
ф = 67.271 0 = 82.287 38 30
ф = 60.42 0 = 54.605 19 34
ф = 66.783 0 = 62.429 27 35
ф = 50.196 0 = 41.611 36 36
Среднее 31±2.514 35.875±1.381
Р 0.207±0.017 0.179±0.007
Таблица 5. Количество выживших дефектов при различных температурах
Направление T = 100 K T = 200 K T = 300 K T = 400 K
ф = 43.314 0 = 70.999 29 26 17 20
ф = 74.023 0 = 71.563 24 26 14 15
ф = 48.309 0 = 59.019 29 19 12 12
ф = 61.289 0 = 30.179 10 14 18 22
ф = 67.271 0 = 82.287 17 16 22 10
ф = 60.42 0 = 54.605 14 8 23 14
ф = 66.783 0 = 62.429 21 9 13 13
ф = 50.196 0 = 41.611 19 23 19 25
Среднее 20.375±2.405 17.625±2.514 17.25±1.436 16.375±1.880
Р 0.203±0.024 0.176±0.025 0.172±0.014 0.164±0.019
Таблица 6. Параметры решетки при различных температурах
Температура Параметры решетки
а A о, A
100 K 3.232 5.170
200 K 3.230 5.173
300 K 3.230 5.176
400 K 3.230 5.178
параметры решетки при различных температурах, которые использовались для построения соответствующих кристаллитов. Параметры решетки были рассчитаны для каждой из рассматриваемых температур методом молекулярной
статики при нулевом давлении. Процент выживших дефектов не превысил 20%.
Стоит отметить, что количество выживших дефектов зависит от параметров решетки кристаллита, т.е. зависит от напряжений в кристал-
ле. Это - отдельный результат, который представляет интерес дальнейшего изучения.
ЗАКЛЮЧЕНИЕ
Кратко сформулируем основные результаты работы.
В работе дан краткий обзор особенностей применения метода молекулярной динамики для моделирования радиационного кристаллических структур.
Путем компьютерного моделирования каскадов атомных смещений методом молекулярной динамики вычислены значения пороговых энергий смещения циркония для различных кристаллографических направлений вылета ПВА.
Получены данные о количестве выживших дефектов при использовании двух разных потенциалов [2, 3]. Результаты хорошо согласуются с проведенными ранее расчетами, однако установлено, что для больших энергий ПВА потенциал [3] дает больший разброс. Количество выживших дефектов при энергиях ПВА 10-20 кэВ не превышало 20%.
Отсутствие зависимости выживших дефектов от температуры подтверждается серией проведенных расчетов в диапазоне температур от 100 до 400 К и хорошо согласуется с данными других работ [4, 10].
В дальнейшем планируется провести расчетные исследования влияния элементов внутренней структуры (границы зерен, выделения 2-х фаз и т.п.) циркония и сплавов на его основе на особенности процесса радиационного повреждения.
Работа выполнена при поддержке Минобрна-
уки в рамках государственного задания на 20122014 гг.
СПИСОК ЛИТЕРАТУРЫ
1. КирсановВ.В. ЭВМ-эксперимент в атомном материаловедении. М.: Энергоатомиздат, 1990. 304 с.
2. Defect, surface and displacement-threshold properties of a-zirconium simulated with many-body potential / G.J. Ackland, S.J. Wooding, D.J. Bacon // Phil. Mag. A, 1995, Vol.71, №3, p.553-565.
3. Mendelev M.I., Ackland GJ. Development of an interatomic for the simulation of phase transformations in zirconium // Phil. Mag. Let., 2007, Vol.87, №5, p.349-359.
4. Molecular dynamics simulations of irradiation cascades in alpha-zirconium under macroscopic strain/ Sali Di, Zhongwen Yao, Mark R. Daymond, Fei Gao // Nuclear Instruments and Methods in Physics Research B 303, 2013, pp. 95-99.
5. Khater H.A., Bacon D.J. // Acta Mater. 58 (2010) 2978.
6. N.d. Diego, A. Serra, D.J. Bacon, Y.N. Osetsky // Model. Simul. Mater. Sci. Eng. 19 (2011) 035003.
7. Finnis M.F., Sinclair J.E. A simple empirical N- Body potential for transition metals // Philos. Mag., A 50, 1984, pp. 45 - 55.
8. Stability of self-interstitial atoms in hcp-Zr / Qing Peng, WeiJi, Hanchen Huang, Suvranu De //Journal of Nuclear Mat, 2012, Vol 429, p.233-236.
9. Charge-optimized many-body (COMB) potential for zirconium / Mark J. Noordhoekm Tao Liang, Zizhe Lu, Tzu-Ray Shan, Susan B. Sinott, Simon R. Phillpot // Journal of Nuclear Mat., 2013, Vol.441, p.274-279.
10. Statistics of primary damage creation in high-energy displacement cascades in copper and zirconium / R.E. Voskoboinikov, Yu.N. Osetsky, D.J. Bacon // Nuclear Instruments and Methods in Physics Research B 242, 2006, pp. 68-70.
MODELING OF THE HCP-ZIRCONIUM BY MOLECULAR DYNAMICS METHOD
© 2013 P.E. Kapustin
Ulyanovsk State University
This paper is devoted to studying the radiation damage of the HCP-Zirconium by molecular dynamics method. We performed calculation of the point defect formation energies and the threshold displacement energy. Simulation of atomic displacement cascades at different temperatures, primary knocked-on atoms energies, and interatomic potentials was carried out. The obtained results are compared with the available literary data. Key words: zirconium, molecular dynamics method, point defects, threshold displacement energy, atomic displacement cascade.
Pavel Kapustin, Graduate Student, Associate Research Fellow of TechnologicalResearchlnstitute. E-mail: [email protected]
