ЖЭТФ, 2023, том 163, вып. 2, стр. 201-213

ОБРАЗОВАНИЕ, ДИФФУЗИЯ И РОСТ ГАЗОНАПОЛНЕННЫХ

ПУЗЫРЬКОВ В γ-УРАНЕ ПРИ ИЗБЫТКЕ МЕЖДОУЗЕЛЬНЫХ

АТОМОВ: СВЯЗЬ МОЛЕКУЛЯРНОЙ ДИНАМИКИ И КИНЕТИКИ

Е. А. Лобашевa,b*, А. С. Антроповa,b**, В. В. Стегайлов^a,b,c

^a Объединенный институт высоких температур Российской академии наук

125412, Москва, Россия

^b Московский физико-технический институт (Национальный исследовательский университет)

141701, Долгопрудный, Московская обл., Россия

^c Национальный исследовательский университет “Высшая школа экономики”

101000, Москва, Россия

Поступила в редакцию 22 августа 2022 г.,

после переработки 14 сентября 2022 г.

Принята к публикации 18 октября 2022 г.

Для описания эволюции ядерных топлив важным процессом является образование газонаполненных на-

норазмерных пузырьков в результате объединения отдельных продуктов деления урана. Теоретическое

описание этого процесса связано со значительными трудностями, так как требует учета в рамках еди-

ной модели как динамики отдельных атомов в решетке, так и кинетики эволюции ансамбля пузырьков.

В данной статье описана попытка построить такую модель, основанную на молекулярно-динамических

(МД) расчетах для пузырьков ксенона в ОЦК-уране в случае избытка междоузельных атомов в кристал-

лической матрице. Анализ основан на МД-моделировании неравновесного процесса образования нанопу-

зырьков ксенона из отдельных атомов Xe, растворенных в кристаллической матрице. Проанализировано

соотношение размера пузырьков и числа атомов газа в них, а также зависимость коэффициента диффу-

зии пузырьков от их радиуса и числа междоузельных атомов в матрице γ-U. Предложена кинетическая

модель эволюции ансамбля пузырьков, позволяющая описать результаты МД-расчетов и экстраполиро-

вать их на большие времена.

DOI: 10.31857/S0044451023020074

газа и нуклеации пузырьков различаются из-за су-

EDN: OQNPWL

щественного различия в размерах атомов газа.

Ключевым этапом выхода газа из материала яв-

ляется его выход на границы зерен. На движение га-

за к границам зерен могут влиять следующие про-

1. ВВЕДЕНИЕ

цессы [3]: диффузия атомов газа в матрице, обра-

Детальное понимание механизмов выхода газо-

зование и рост пузырьков за счет поглощения ва-

образных продуктов деления (преимущественно ксе-

кансий или испускания междоузельных атомов, по-

нона) из ядерных топлив для математического моде-

глощение газа пузырьками и перерастворение газа,

лирования этого процесса представляет собой важ-

а также два механизма слияния пузырьков — ми-

ную современную научно-техническую проблему [1,

грация и коалесценция и созревание по Оствальду

2]. С качественной точки зрения эти процессы ана-

(последний механизм играет роль при высоких тем-

логичны процессам транспорта гелия в металлах,

пературах и низких концентрациях газа).

которые изучаются для конструкционных матери-

Многие вычислительные подходы учитывают

алов ядерных и термоядерных реакторов [3]. При

только диффузию атомов газа в матрице, захват

этом, конечно, конкретные механизмы диффузии

(trapping) атомов газа неподвижными пузырьками

и перерастворение (resolution) атомов газа из пу-

зырьков в матрицу в результате воздействия тре-

^* E-mail: lobashev.ea@phystech.edu

^** E-mail: antropov@phystech.edu

ков осколков деления [4,5]. Игнорирование иных ме-

201

Е. А. Лобашев, А. С. Антропов, В. В. Стегайлов

ЖЭТФ, том 163, вып. 2, 2023

ханизмов выхода газа из объема зерен материала

сий. Однако результаты Верма с соавторами [14] по

топлива обусловлено, вообще говоря, как нехваткой

расчету скорости выхода ксенона из диоксида ура-

экспериментальных данных, так и сложностью тео-

на в рамках многомасштабной модели, учитываю-

ретического описания соответствующих процессов.

щей ускоренное движение пузырьков в градиенте

Так, например, зачастую используется упрощенное

концентрации вакансий [15], показывают очень зна-

описание, согласно которому в результате перерас-

чительные расхождения с результатами специально

творения появляются пузырьки одинакового разме-

поставленных экспериментов.

ра, хотя имеются указания на то, что учет распреде-

С помощью разработанного ранее метода нерав-

ления пузырьков по размерам может повысить точ-

новесной молекулярной динамики (МД) для моде-

ность подобных моделей [6].

лирования диффузии пузырьков [16,17] недавно на-

Некоторые топливные коды (например, код

ми были получены результаты [18] для диоксида

MFPR, разрабатывающийся в ИБРАЭ РАН [7])

урана, свидетельствующие о существовании меха-

явно включают в математическую модель мигра-

низмов сверхбыстрой диффузии газонаполненных

цию и коалесценцию пузырьков. В этом случае

нанопузырьков, которые пока не учитываются в

отдельного внимания заслуживают коэффициенты

топливных кодах.

диффузии газонаполненных пузырьков.

Как было указано в работе [19], эволюция попу-

Для отдельных задач, связанных с поведением

ляции пузырьков через миграцию и коалесценцию

газовых пузырьков, были предложены аналитиче-

подчиняется уравнениям, схожим с таковыми для

ские решения. В работах Рязанова, Волкова и Вос-

коагуляции коллоидных частиц и может быть опи-

кобойникова была развита общая теория образова-

сана моделью Чандрасекара [20].

ния пузырьков в перенасыщенных растворах вакан-

МД-моделирование активно используется для

сий, междоузельных атомов и атомов газа [8-10]. В

исследования поведения пузырьков Xe в UO₂. Мур

частности, было показано, что высокая плотность

с соавторами в работе [21] установили, что кластер

газа в пузырьках заметно снижает размер крити-

атомов Хе формируется из отдельных атомов Хе, за-

ческого зародыша и барьер его образования [8]. На

нимающих вакансии Шоттки, чему способствует пе-

основе этого результата в работе Реста [11] был про-

ресыщение вакансий Шоттки в UO₂. Мерфи с соав-

анализирован механизм нуклеации газовых пузырь-

торами в работе [22] изучали свободную энергию Хе

ков одновременно из нескольких атомов газа (в про-

в точечных дефектах и пузырьках в UO₂. Джелеа

тивоположность общепризнанному ранее двухатом-

с соавторами [23] также изучали поведение пузырь-

ному механизму). В работе Нуаро показана связь

ков Хе в UO₂ и их влияние на распухание топли-

равновесных концентраций газа в пузырьке с его

ва. Лиу с соавторами [24] моделировали нуклеацию

размером и концентрациями дефектов в объеме мат-

пузырьков ксенона и гелия в UO₂. Кинетика роста

рицы [12]. Вещуновым и Шестаком была развита

ксеноновых пузырьков в диоксиде урана и испуска-

детальная теория диффузии газонаполненных пу-

ние ими дислокационных петель изучалось методом

зырьков в диоксиде урана [13].

МД в [25]. Другим методом моделирования разви-

Важной подзадачей является изучение транспор-

тия газовой пористости в материалах является ме-

та газа в объеме матрицы без радиационного сти-

тод фазового поля [26], к которому, однако, нужно

мулирования (т. е. при отсутствии в материале тре-

относится как к качественному методу.

ков осколков деления атомов урана). В этом слу-

Концентрации Xe в пузырьках радиусом

чае радиационное перерастворение газа перестает

500 нм были экспериментально измерены [27] и ока-

играть роль и определяющей становится диффузия

зались выше равновесных значений, что можно объ-

пузырьков. В такой постановке результаты теоре-

яснить механизмом испускания дислокационных пе-

тических расчетов можно сравнить с результатами

тель при росте пузырьков [28]. Поэтому отдельный

экспериментов по выходу газа из облученного мате-

интерес представляет рост пузырьков в условиях из-

риала при отжиге, как, например, в [14]. И в этом

быточной концентрации междоузельных атомов.

случае вопрос о скорости выхода газа из топлива

Обсуждаемые механизмы гораздо меньше иссле-

также нельзя считать решенным. В работе Эван-

дованы в других видах топлива. Металлическое топ-

са [15] указывается на то, что выход газа из объема

ливо UMo — один из перспективных кандидатов на

зерна происходит гораздо быстрее, чем это предска-

роль топлива для реакторов на быстрых нейтронах.

зывают теоретические модели. Там же предполага-

Пузырьки ксенона в сплаве UMo уже исследовались

ется, что ускоренный выход газа связан с движени-

методами атомистического моделирования [29, 30],

ем пузырьков в поле градиента концентрации вакан-

однако в них не рассматривались начальные этапы

202

ЖЭТФ, том 163, вып. 2, 2023

формирования пузырьков. Экспериментальная ра-

печивают зачастую более высокую точность, их вы-

бота [31] по наблюдению пузырьков, образующихся

числительная сложность на 1-2 порядка выше.

при бомбардировке монокристалла молибдена иона-

Начальное состояние системы — ОЦК-решетка

ми ксенона дала возможность верифицировать мо-

из атомов урана, представляющая собой куб со сто-

дели диффузии ксенона в молибдене [32] (в данной

роной 20 постоянных решетки. Из решетки уда-

модели ключевая роль отводилась радиационно сти-

лялись атомы урана (с номерами, выбранными по

мулированной диффузии). В работе [33] моделиро-

случайному принципу), чтобы обеспечить некото-

валась нуклеация пузырьков ксенона в молибдене

рую начальную концентрацию вакансий, а также

из растворенных атомов газа, находящихся в вакан-

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

сии. В этой работе было показано наличие порого-

нона (после этого проводилась минимизация энер-

вой концентрации ксенона, необходимой для образо-

гии системы). Наиболее интересные результаты уда-

вания пузырьков, а также механизм ускорения диф-

лось получить для 127 атомов ксенона. Система при-

фузии вакансий вблизи пузырька.

водилась к состоянию равновесия в течение 1 нс

В данной работе рассматривается чистый ОЦК-

при помощи термостата Нозе-Гувера при темпера-

уран как прототип топлива UMo. Особенностью

туре 1400 K и нулевом давлении. Такая темпера-

ОЦК-урана является аномально низкая энергия об-

тура была выбрана из методических соображений:

разования междоузельных атомов [34] и их пре-

при больших температурах происходит плавление,

имущественная роль в самодиффузии. Это приво-

а при меньших диффузия происходит медленнее.

дит, например, к необычному механизму диффу-

Несмотря на то, что некоторые атомы урана бы-

зии нанополостей [17]. Методами классической мо-

ли удалены из решетки, в системе сразу в начале

лекулярной динамики моделируется эволюция газо-

МД-расчета появлялись междоузельные атомы. Это

вых пузырьков в решетке ОЦК-урана из растворен-

связано с тем, что атомы ксенона вытесняли атомы

урана в междоузельные положения и занимали их

ных в ней атомов ксенона. Обсуждаются как тер-

модинамические характеристики равновесной попу-

место за время порядка 10 пс (энергия образования

междоузельных атомов U меньше, чем междоузель-

ляции пузырьков, так и кинетика их объединения

и роста. В разд. 2 описываются условия моделиро-

ных атомов Xe, поэтому междоузельные атомы Xe в

вания. В разд. 3 рассматривается связь между раз-

расчетах не наблюдались). За это же время группы

мерами пузырьков, количеством атомов ксенона в

из 2-7 атомов ксенона, начальные положения кото-

них и концентрацией точечных дефектов. В разд. 4

рых оказались рядом, сформировали пузырьки.

вычисляются коэффициенты диффузии пузырьков

Часть появившихся междоузельных атомов ан-

небольших размеров. В разд. 5 рассматривается тео-

нигилировала со свободными вакансиями, поэтому

ретическая модель эволюции популяции газовых пу-

количество междоузельных атомов в начальный мо-

зырьков, основанная на подходе Чандрасекара [20]

мент времени примерно равно количеству атомов

и предсказания этой модели сравниваются с резуль-

ксенона минус начальное количество вакансий. Наи-

татами моделирования.

более подробно был изучен случай с 67 начальными

вакансиями, потому что в этом случае концентрация

междоузельных атомов после релаксации оказыва-

ется близка к равновесной при данной температуре,

2. МОДЕЛЬ И МЕТОД РАСЧЕТА

а именно 127 - 67 = 60 атомов на ячейку.

Также были проведены расчеты, в которых на-

В работе проводилось моделирование эволюции

чальное количество вакансий было выше (до 257). В

решетки урана с растворенными в ней в начальный

этом случае даже после того, как атомы ксенона вы-

момент атомами ксенона. Был использован EAM-

теснили атомы урана, свободных вакансий осталось

потенциал для системы U-Mo-Xe [35]. Этот потен-

больше, чем междоузельных атомов.

циал был разработан для моделирования процес-

Вакансии и междоузельные атомы определя-

сов с участием ксенона в сплаве U-Mo и уже про-

лись в пакете OVITO [42] (рассматривались ячей-

шел апробацию в исследованиях различных науч-

ки Вигнера-Зейтца, окружающие идеальные поло-

ных групп (см., например, [30, 36-39]). Модель хо-

жения узлов решетки, вакансией считалась ячейка,

рошо описывает и свойства чистых U, Mo и Xe.

в которой не было атомов, а междоузельным ато-

Хотя активно развивающиеся в настоящее время

мом — ячейка, в которой находились два атома).

машинно-обучаемые и нейросетевые модели потен-

Далее система рассчитывалась в NVE-ансамбле

циалов для многокомпонентных систем [40,41] обес-

с шагом 1 фс. Постоянная решетки a

= 3.555Å

203

Е. А. Лобашев, А. С. Антропов, В. В. Стегайлов

ЖЭТФ, том 163, вып. 2, 2023

12,5

7,5

2,5

Csia = (5.6÷8.3)⋅ 10^-4Å^-3

Csia = (1.4÷2.8)⋅ 10^-4Å^-3

100

атомы ксенона

Рис. 2. Зависимость радиуса пузырька от количества ато-

мов ксенона n для двух концентраций междоузельных ато-

мов урана. Линия соответствует формулам (2) и (1

Была рассмотрена гипотеза, что размер пузырь-

ка для фиксированного n (n — число атомов Xe в

пузырьке) может зависеть от концентрации меж-

доузельных атомов C_sia (self-interstitial atoms), ко-

торая существенно меняется со временем. Поэто-

му сначала рассматривались только пузырьки, су-

Рис. 1. Моделируемая система после 30 нс МД-расчета.

ществовавшие в начале моделирования, а потом

На верхней картинке показаны все атомы (красные — U,

только существовавшие в конце (два набора дан-

желтые — Xe), а на нижней — результат анализа той же

ных на рис. 2). Рассмотренные в работе значения

структуры (атомы урана не показаны). Голубым цветом

концентраций междоузельных атомов C_sia соответ-

обозначены пустые вакансии, желтым — вакансии, внутри

ствуют отношениям концентрации C_sia к концентра-

которых содержится атом ксенона

ции атомов урана (равной 440 · 10-4 Å-3) порядка

0.003-0.02. По полученному графику можно сделать

соответствует нулевому давлению при температу-

вывод, что зависимость размера пузырьков от кон-

ре 1400 K. Моделирование проводилось в течение

центрации газа практически не меняется при изме-

100 нс. Для системы с 67 вакансиями было рас-

нении концентрации междоузельных атомов урана.

считано 20 траекторий по 100 нс каждая (траек-

тории различаются только направлениями скоро-

Для пузырьков размера n = 5 и более атомов ксе-

стей атомов в начальный момент времени, что до-

нона количество вакансий, приходящихся на один

статочно для стохастизации, приводящей к их ста-

атом ксенона, практически не зависит от размера

тистической независимости [43]). Для анализа ме-

и равно примерно 2.7 (здесь мы называем объемом

тодом Вигнера-Зейтца было необходимо избавить-

вакансии объем одной ячейки Вигнера-Зейца ОЦК-

ся от слишком больших тепловых колебаний, по-

фазы урана). А для меньших пузырьков это отно-

этому положения атомов усреднялись по 1000 МД-

шение меньше (рис. 3), например, пузырек с одним

шагам (рис. 1). В рассматриваемой МД-модели за

атомом ксенона состоит из одной вакансии. Таким

100 нс давление и температура в системе меняют-

образом, зависимость количества вакансий, состав-

ся незначительно (давление возрастает примерно до

ляющих пузырек, от количества атомов ксенона мо-

6000 бар, а температура поднимается на 80 К). Рас-

жет быть приближенно описана эмпирической фор-

четы проводились с помощью программного ком-

мулой

плекса LAMMPS [44] на суперкомпьютере “Десмос”

в ОИВТ РАН [45].

N_vac/n = (2.7 - 2.1 exp(-0.27n)).

(1)

204

ЖЭТФ, том 163, вып. 2, 2023

3. СВОЙСТВА ПОПУЛЯЦИИ ПУЗЫРЬКОВ

3.1. Соотношение размера пузырьков и

числа атомов газа

2.7

Сразу после термостатирования большинство

атомов ксенона находилось в виде атомов замеще-

ния в решетке; после этого они начали объединять-

ся в пузырьки. Процессы перехода атомов Xe между

пузырьками (созревание по Оствальду) в данной си-

стеме практически не наблюдались.

На рис. 4 показано, что число междоузельных

атомов урана растет с течением времени в результа-

те укрупнения газовых пузырьков.

Для всех пузырьков, содержащих одинаковое ко-

Csia = (1.4÷8.3)⋅ 10^-4Å^-3

личество атомов ксенона, было найдено среднее ко-

100

личество вакансий N_vac, составляющих их. На рис. 2

атомы ксенона

показана зависимость эффективного радиуса пу-

Рис. 3. Количество вакансий, приходящихся на один атом

зырька от количества атомов ксенона. За радиус

ксенона, для пузырьков различного размера. Линия соот-

был принят радиус шара, объем которого равен сум-

ветствует формуле (1)

марному объему вакансий:

)1/3

⁽3N_vac

300

R=a

(2)

8π

Отдельно стоит сказать про расчеты с повышен-

250

ной начальной концентрацией вакансий. Выше бы-

ло сказано, что даже после того, как атомы ксенона

200

вытесняют атомы урана в междоузельные положе-

Nvac(0) = 67

ния, в системе остаются свободные вакансии. Оказа-

Nvac(0) = 257

лось, что в таких расчетах все свободные вакансии

150

почти сразу объединялись с пузырьками, и в систе-

ме все равно появлялось некоторое количество меж-

100

доузельных атомов. При этом количество вакансий

на атом ксенона в пузырьке было больше, чем по-

лучается по формуле (1). Но через некоторое время

в результате объединения пузырьков количество ва-

кансий на атом ксенона достигло такого же значе-

ния (т. е. в некотором смысле произошло “насыще-

100

ние” пузырьков), и после этого концентрация меж-

время, нс

доузельных атомов начала расти, как в предыдущем

Рис. 4. Зависимость количества междоузельных атомов

случае (рис. 4).

от времени. Кривыми показаны зависимости, получен-

ные в результате решения кинетической модели (уравне-

3.2. Давление газа в пузырьках

ния (20) и (22))

В континуальном пределе зависимость равновес-

ного давления газа в пузырьке от его радиуса опре-

В результате того, что в больших пузырьках на

деляется формулой Лапласа

один атом ксенона приходится больше вакансий,

2γ

при объединении пузырьков некоторые атомы ура-

P =

(3)

на вытесняются в междоузельные положения. Этим

и объясняется зависимость концентрации междо-

где γ — коэффициент поверхностного натяжения.

узельных атомов от времени, показанная на рис. 4.

В связи с тем, что пузырьки находятся в кристал-

лической матрице, следует также учесть поправку,

205

Е. А. Лобашев, А. С. Антропов, В. В. Стегайлов

ЖЭТФ, том 163, вып. 2, 2023

150

связанную с выдавливанием газом междоузельных

L = 0.5 Å

атомов (или дислокационных петель, представляю-

L = 0

щих собой кластеры из нескольких междоузельных

125

P= 2γ/R

атомов, имеющих форму плоских “дисков”) [46]:

2γ

E_i/l

100

P =

(4)

V_i/l

где Ei/l — энергия образования междоузельного

атома или дислокационной петли, а V_i/l — соответ-

ствующий объем, освобождающийся при их испус-

кании.

Для твердого вещества нет простого способа

определения коэффициента поверхностного натяже-

ния, но для количественного анализа нам будет до-

статочно оценки. Согласно [47] для ОЦК- и ГЦК-

решеток при нулевой температуре есть универсаль-

R, Å

ная связь между энергией атома в решетке E_coh и

коэффициентом поверхностного натяжения (выра-

Рис. 5. Зависимость давления ксенона от его концентра-

женного в эВ/атом)

ции. Линия определяется формулой (3), точки соответ-

ствуют данным МД и соотношениям (8) и (9)

γ(0) = 0.134E_coh.

(5)

Для ОЦК-урана в рамках используемой мо-

Можно также учесть, что эффективный ради-

дели потенциала EAM было рассчитано, что

ус пузырька, использующийся для расчета концен-

E_coh = 4.147 эВ/атом.

трации газа в нем, может отличаться от опреде-

В [48] предлагается следующая зависимость по-

ления (2). Из-за того, что среднее расстояние U-

верхностного натяжения от температуры:

Xe несколько больше, чем среднее расстояние U-U,

целесообразно ввести поправку на объем пузырька

γ(T ) = γ(0) - T S,

(6)

для расчета концентрации газа в нем. Толщина соот-

ветствующего “исключенного слоя” может быть оце-

где принимается энтропийный вклад S = 1.0k_B.

нена как L = 0.5Å при помощи сравнения функций

Согласно формулам (5) и (6) γ = 0.049 эВ/Å2

радиального распределения U-U и U-Xe. Тогда в

(переход от эВ/атом к эВ/Å2 делается с использо-

формуле (8) R следует заменить на (R - L):

ванием средней площади, приходящейся на атом по-

верхности пузырька). Кривая, описываемая форму-

лой (3), показана на рис. 5.

C =

(9)

4π(R - L)³

Для того чтобы сравнить результаты молекуляр-

ной динамики с формулами (3) и (4), необходимо

Точки, рассчитанные с этой поправкой, также изоб-

знать зависимость давления ксенона от концентра-

ражены на рис. 5.

ции (плотности). Для этого была проведена серия

Видно, что серые точки (L = 0) на рис. 5 выхо-

МД-расчетов ячейки с чистым ксеноном (2000 ато-

дят на константу P = 25 кбар. Можно ожидать, что

мов). Получено, что

черные точки (L = 0.5Å) выйдут на ту же констан-

ту, только при больших радиусах R ≫ L. Если пред-

P =a_XeCbXe,

(7)

положить механизм роста пузырьков за счет испус-

кания одиночных междоузлий, то асимптотическое

где C

— концентрация ксенона в

Å-3,

давление должно составлять E_i/(a³/2) ∼ 50 кбар,

a_Xe

2.54 · 10¹⁰ бар, b_Xe

= 3.38. Эта оцен-

где E_i ∼ 0.7 эВ [34], a³/2 — атомный объем U в ОЦК-

ка достаточно точно выполняется в диапазоне

матрице. Таким образом, асимптотическое давление

C = (0.015-0.033) Å-3. Концентрация определялась

газа в пузырьках в МД-модели оказывается пример-

как

но в два раза ниже, чем если бы пузырек расши-

C =

(8)

4πR³

рялся за счет одиночных междоузлий. При испус-

На рис. 5 показаны точки, полученные по формулам

кании дислокационных петель энергия образования

(7) и (8).

E_l дислокационной петли из N_l междоузлий может

206

ЖЭТФ, том 163, вып. 2, 2023

быть ниже, чем N_lE_i, за счет изменения энергии

n=1, поделено на 5

междоузлий при объединении их в дислокационную

n=6

петлю. Однако оценить эту энергию намного слож-

n=10

n=25

нее, чем энергию образования междоузельного ато-

1,5

ма E_i.

4. ДИФФУЗИЯ ПУЗЫРЬКОВ

Для параметризации кинетической модели, опи-

сывающей рост пузырьков (см. разд. 5), необходи-

мо определить коэффициенты диффузии пузырьков

различного размера. Было принято, что коэффици-

0,5

ент диффузии пропорционален концентрации меж-

доузельных атомов [18]. Для проверки этого пред-

положения была отдельно рассчитана зависимость

коэффициента диффузии пузырька из 10 атомов

от концентрации междоузельных атомов. Ячейка со

2×10^-4

4×10^-4

6×10^-4

8×10^-4

стороной 10 постоянных решетки с единственным

Сsia, Å-3

пузырьком моделировалась в течение 260 нс. На-

Рис. 6. Зависимость коэффициента диффузии пузырьков

чальная концентрация междоузлий варьировалась

различных размеров от концентрации междоузельных ато-

в расчетах, в процессе моделирования она точно

мов урана

определялась методом Вигнера-Зейтца и усредня-

лась за все время моделирования. Остальные пара-

метры были такими же, как в основных 20 расчетах.

в ячейке с множеством пузырьков

На рис. 6 показаны полученные зависимости (для

8×10⁵

в ячейке с единственным пузырьком

ясности отметим, что 1Å2/нс = 10-7 см²/с). Пред-

положение о пропорциональной зависимости оказы-

вается разумным, хотя для некоторых размеров на-

6×10⁵

блюдаются существенные отклонения от нее.

Таким образом, мы можем вычислять не сам ко-

эффициент диффузии, а отношение D/C_sia .

Расчет коэффициента диффузии газовых пу-

зырьков проводился по тем же 20 траекториям сле-

дующим образом. Для каждого пузырька было най-

дено смещение за 1 нс (если за это время он не объ-

единился с другим пузырьком), коэффициент диф-

2×10⁵

фузии был определен по формуле

〈Δr²〉

(10)

6Δt

и нормирован на среднюю концентрацию междо-

время, нс

узельных атомов на протяжении данной одной на-

носекунды. На рис. 7 показана зависимость среднего

Рис. 7. Средний квадрат смещения пузырька из 10 ато-

квадрата смещения пузырька от времени.

мов, деленный на концентрацию междоузельных атомов,

Получившиеся значения коэффициента диффу-

в зависимости от времени. Для сравнения приведен также

зии показаны на рис. 8. Минимумы около n = 7 и

график, полученный в результате моделирования ячейки

с одним пузырьком

n = 25 атомов объясняются правильной геометриче-

ской формой соответствующих пузырьков (рис. 9).

Для пузырьков из 30 и более атомов Xe коэффици-

что для больших пузырьков D ∝ R-3 ∝ n-1 или

енты диффузии были рассчитаны с меньшей точ-

D ∝ R-4 ∝ n-3 в зависимости от механизма диф-

ностью, поэтому для них они были приближены

фузии (объемный или поверхностный соответствен-

асимптотической зависимостью. Можно ожидать,

но [13]). По полученным МД-данным невозможно

207

Е. А. Лобашев, А. С. Антропов, В. В. Стегайлов

ЖЭТФ, том 163, вып. 2, 2023

R, Å

для обоих случаев. Для поверхностного механизма

02 3

асимптотика была приближена следующей форму-

лой:

)

190000

⁽A

n ≥ 30,

(11)

C_sia

n4/3

нс

а для объемного — следующей:

49000

⁽A⁵)

n ≥ 30.

(12)

C_sia

n нс

5. КИНЕТИЧЕСКАЯ МОДЕЛЬ

5.1. Формулировка модели

В данном разделе описывается модель диффу-

зии и нуклеации пузырьков, основанная на подхо-

де Чандрасекара [20]. В работе [20] рассматривает-

ся модель для коллоидных частиц в растворе, мы

кол-во атомов ксенона

же используем аналогичный математический под-

Рис. 8. Зависимость коэффициента диффузии от разме-

ход для газовых пузырьков.

ра при концентрации междоузельных атомов 6.9·10-4Å-3.

Будем считать, что пузырьки не взаимодейству-

Коэффициент диффузии при любой другой концентрации

ют, находясь на расстояниях больше R, а сблизив-

можно получить, учитывая, что D ∝ Csia

шись на это расстояние, они мгновенно объединяют-

ся в один.

Сначала будем рассматривать один неподвиж-

ный пузырек размера n, находящийся в начале ко-

ординат, вокруг которого равномерно распределены

пузырьки одинакового размера m с начальной кон-

центрацией C_m. Если решить уравнение диффузии

= DΔC

(13)

с граничными условиями

{

C =C

|r| > R, t = 0,

(14)

C = 0,

|r| = R, t > 0,

то получится следующее:

(

))

⁽r-R

C =C_m

r - R + R·erf

√

(15)

Отсюда частота объединения такого пузырька с

Рис. 9. Сверху: пузырьки из 15 вакансий (7 атомов ксено-

остальными будет равна

на) и 62 вакансий (25 атомов ксенона). Снизу: идеальный



(

)

dC

пузырек из 65 вакансий и изображение ромбического до-



J = 4πR²D

= 4πDRC_m 1 +

√

(16)



декаэдра для сравнения

πDt

r=R

где вторым слагаемым обычно можно пренебречь,

однозначно определить, какой из этих механизмов

если рассматривается достаточно большой отрезок

преобладает в данном случае, так как они пример-

времени:

но с одинаковой точностью описывают зависимость

D(n) при n ≥ 30. Поэтому были проведены расчеты

t≫

(17)

R²

208

ЖЭТФ, том 163, вып. 2, 2023

Если учесть, что рассматриваемый пузырек тоже

пузырьки, а также начальной концентрацией вакан-

может двигаться, то в формуле (16) вместо D бу-

сий.

дет коэффициент взаимной диффузии

Результаты решения кинетических уравнений

и результаты моделирования показаны на рис. 10.

D_nm = D_n + D_m.

(18)

Для определения количества и размеров пузырьков

в моделировании за пузырек считался кластер ато-

Теперь можно записать частоту объединений пу-

мов ксенона, т. е. группа атомов, находящихся друг

зырьков двух размеров:

от друга на расстояниях не более 4˚. Концентрация

J_nm = 4π(D_n + D_m)R_nmC_nC_m

(19)

пузырьков находилась как количество пузырьков в

ячейке, деленное на ее объем, т. е. на 3.6 · 10⁵ Å³.

и искомую систему кинетический уравнений:

Сравнение результатов МД-расчетов (усредненных

∑

по 20 траекториям) и результатов кинетической мо-

dC_n

= -4πC_n C_mR_nmD_nm +

дели показывает, что кинетическая модель обеспе-

(20)

чивает адекватное описание МД-данных.

∑

+ 4π

C_mC_kR_mkD_mk.

Преимущество кинетической модели в том, что

m+k=n

с ее помощью можно рассчитывать распределение

газовых пузырьков по размерам на временах, су-

Здесь первое слагаемое означает исчезновение пу-

щественно превышающих времена, доступные пря-

зырьков размера n за счет объединения с другими,

мым МД-расчетам. В качестве примера был прове-

а второе — появление пузырьков размера n за счет

ден расчет, предсказывающий распределение по раз-

объединения меньших пузырьков, сумма размеров

мерам через одну микросекунду после начала рас-

которых равна n.

чета (рис. 11). Поскольку в системе при этом пре-

На основе наших расчетов можно проверить точ-

обладают пузырьки больших размеров, для кото-

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

рых коэффициенты диффузии определяются с по-

тимальные параметры в модели.

мощью асимптотической аппроксимации, были рас-

смотрены оба варианта аппроксимации (D ∝ R-3

5.2. Параметризация, проверка и

и D ∝ R-4). Рис. 11 показывает, что на масштабе

применение модели

1 мкс зависимость от выбора асимптотики невели-

В качестве радиуса объединения было выбрано

ка (влияние заметнее на больших временах). Мож-

но заметить, что распределение, предсказанное ки-

R_nm = R_n + R_m + b,

(21)

нетической моделью для времени 1 мкс, имеет пики

около пузырьков с низким коэффициентом диффу-

где b = 6Å было подобрано таким образом, чтобы

зии.

средний размер пузырька через 20 нс после нача-

На рис. 12 показано, как менялась концентрация

ла расчета совпадал с результатами моделирования.

пузырьков некоторых размеров в течение 1000 нс.

Радиус R_n определялся по формулам (2) и (1).

Видно, что имеет место максимум концентрации для

Для решения уравнений (20) был использован

пузырьков каждого размера.

метод Рунге-Кутты 4-го порядка. Учитывались пу-

зырьки размером до 240 атомов. Было принято, что

в начальный момент весь ксенон находится в пу-

6. ОБСУЖДЕНИЕ

зырьках из одного атома (т. е. отдельные атомы Xe

в позициях замещения в узлах ОЦК-решетки).

Рассмотренный в работе относительно неболь-

В кинетические уравнения входит коэффициент

шой размер МД-модели был ориентирован на повы-

диффузии, а он зависит от концентрации междо-

шение физического времени в МД-расчетах. Высо-

узельных атомов. Она рассчитывалась как

кая начальная концентрация растворенного в ОЦК-

∑

4π

уране ксенона позволила ускорить процесс образо-

C_sia = -C_vac(0) +

R3i,

(22)

a³

вания пузырьков и уменьшить времена ожидания

их последующего взаимодействия, что облегчает по-

где R_i — радиус i-го пузырька, a³/2 — объем, при-

лучение кинетической информации из МД [49]. Во-

ходящийся в решетке на один атом. Таким обра-

обще говоря, в дальнейшем планируется рассмот-

зом, количество междоузельных атомов определяет-

реть вопрос влияния начальной концентрации га-

ся суммарным количеством вакансий, образующих

за. Современные суперкомпьютерные системы да-

209

ЖЭТФ, вып. 2

Е. А. Лобашев, А. С. Антропов, В. В. Стегайлов

ЖЭТФ, том 163, вып. 2, 2023

R, Å

04 67

03 4

0,05

0,4

решение кинетических уравнений

атомистическое моделирование

0,3

0,04

0,2

0,1

0,03

10 нс

100

150

200

0,02

R^-3

R^-4

12,5

0,01

мкс

7,5

100

150

200

250

20 нс

размер пузырька (кол-во атомов ксенона)

2,5

Рис. 11. Распределение пузырьков по размерам через 1

мкс после начала расчета, полученное решением кинети-

ческих уравнений

кол-во атомов Xe:

нс

7 (концентрации поделены на 5)

0,8

0,6

1,5

1,25

0,4

50 нс

0,75

0,2

0,5

0,25

200

400

600

800

1000

время, нс

кол-во атомов ксенона в пузырьке

Рис. 12. Изменение концентрации пузырьков с 1, 7, 25 и

Рис. 10. Распределение пузырьков по размерам в различ-

30 атомами Xe в течение 1000 нс

ные моменты времени для 20 траекторий

для γ-U характерна аномально низкая энергия об-

разования междоузельных атомов E_i ∼ 0.7 эВ и, со-

ют возможность моделировать МД-модели из де-

ответственно, аномально высокая равновесная кон-

сятков миллионов атомов с аналогичными EAM-

центрация междоузельных атомов [34]. Следует от-

потенциалами на сотни наносекунд [50], особенно с

метить, что эта особенность γ-U не позволяет при-

использованием графических ускорителей [51, 52].

менять для описания газовых пузырьков в γ-U об-

На первый взгляд, рассмотренная в данной ра-

щую теорию [8-10], разработанную в приближении

боте ситуация избытка междоузельных атомов яв-

(справедливом для большинства металлов), что рав-

ляется достаточно специальным случаем неравно-

новесная концентрация междоузельных атомов на-

весной концентрации точечных дефектов, который

много меньше, чем концентрация вакансий.

может реализовываться только в особенных усло-

Детальный анализ давления газа в пузырь-

виях (например, под облучением). Однако именно

ках показал, что на рассмотренных в работе МД-

210

ЖЭТФ, том 163, вып. 2, 2023

траекториях в газовых пузырьках всех размеров,

облегчает перестройку поверхности пузырьков, что

начиная с наименьших, давление газа выше равно-

обеспечивает колоссальное (на порядки) ускорение

весного лапласовского давления. Хотя существова-

скорости диффузии пузырьков в кристаллической

ние пересжатых (over-pressurized) газовых пузырь-

матрице. Нами установлен факт пропорционально-

ков в топливе U-Mo описано в литературе (см., на-

сти D ∝ C_sia для большинства размеров газовых

пример, [53]), информации о механизме формиро-

пузырьков.

вания таких пузырьков и давлении газа в них все

В то же время результаты работы показывают,

еще недостаточно. Интересным результатом данной

что маленькие пузырьки с определенным числом

работы является существование асимптотики соот-

атомов газа (n = 7 и n = 25) демонстрируют ано-

ношения числа вакансий, составляющих пузырек,

мально низкую скорость диффузии. Причина это-

к числу атомов газа в нем, N_vac/n, для пузырьков

го состоит в особенной регулярной форме подоб-

небольшого размера (R < 1 нм), начиная с n ∼ 20.

ных пузырьков, соответствующей правильным мно-

В недавней работе Билера с соавторами [37] были

гогранникам, в результате чего сильно снижается

рассмотрены более крупные пузырьки ксенона раз-

число подвижных поверхностных атомов и, следова-

мером 3-8.5 нм в сплаве U-Mo (в рамках точно та-

тельно, эффективно снижается и поверхностная са-

кой же модели межатомного потенциала [35], что

модиффузия, обеспечивающая диффузию пузырь-

используется в этой работе), и авторы делают за-

ков. Некоторые причины несферичности полостей в

ключение о возрастании соотношения N_vac/n для

ОЦК-металлах были недавно проанализированы в

пузырьков этих размеров. Асимптотическое значе-

работе Назарова с соавторами [55]. Низкая подвиж-

ние давления 25 кбар ∼ E_i/a³ в пузырьках, рассмот-

ность пузырьков указанных размеров обеспечивает

ренных нами в данной работе, может быть объясне-

их высокую концентрацию на больших временах,

но механизмом выдавливания давлением газа ком-

что демонстрирует решение кинетической модели,

плексов междоузельных атомов урана с поверхности

параметризованной по МД-данным. Наличие высо-

пузырька в объем ОЦК-матрицы. Детали подобно-

кой концентрации подобных малоподвижных сверх-

го механизма требуют специального исследования,

малых пузырьков может быть важно с точки зре-

аналогичного, например, работе по МД-анализу ис-

ния интерпретации экспериментальных данных и

пускания дислокационных петель содержащими ге-

построения моделей выхода газа из топлива. Отме-

лий нанометровыми пузырьками в вольфраме [54].

тим, что наличие большого числа подобных пузырь-

Для того чтобы применить разработанную мо-

ков может ускользать из рассмотрения, поскольку

дель к системе с избыточной концентрацией вакан-

предел разрешения основного экспериментального

сий, необходимо решить ряд вопросов. Во-первых,

метода анализа пузырьков — просвечивающей элек-

предположение о том, что коэффициент диффузии

тронной микроскопии (TEM) — соответствует 1 нм.

пропорционален концентрации междоузельных ато-

мов, не позволяет проводить расчеты при нулевой

концентрации междоузельных атомов. Во-вторых,

7. ЗАКЛЮЧЕНИЕ

при низких концентрациях междоузельных атомов

важно учитывать зависимость соотношения N_vac/n

Проанализирована модель образования газовых

не только от n, но и от C_sia . На рис. 4 для случая с

нанопузырьков из ксенона, растворенного в матри-

257 начальными вакансиями решение было построе-

це ОЦК-урана при избытке междоузельных атомов.

но, только начиная с момента 30 нс, так как до этого

Показано, что для пузырьков размером R ∼ 5-10Å

момента концентрация междоузельных атомов была

= 20-100 атомов Xe) характерна одинаковая

слишком мала. Кроме того, для описания случая с

концентрация газа. Соответствующее максимальное

пониженной концентрацией междоузельных атомов

предельное давление газа заметно выше равновес-

соотношение (1) потребовало уточнения.

ного лапласовского и определяется механизмом вы-

давливания газом атомов урана с поверхности пу-

Результаты данной работы свидетельствуют об

зырька в объем ОЦК-матрицы.

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

Получены коэффициенты диффузии пузырьков

фузии пересжатых газонаполненных пузырьков, что

в рассмотренном диапазоне размеров. Оказалось,

говорит об универсальности подобного механизма,

что зависимость коэффициента диффузии от разме-

уже обнаруженного нами недавно в диоксиде ура-

ра для пузырьков из 1-30 атомов имеет существен-

на [18]. Наличие междоузельных атомов в припо-

но немонотонный вид. Обнаружено существование

верхностном слое газовых пузырьков существенно

211

Е. А. Лобашев, А. С. Антропов, В. В. Стегайлов

ЖЭТФ, том 163, вып. 2, 2023

малоподвижных пузырьков, форма которых соот-

11.

J. Rest, J. Nucl. Materials 402, 179 (2010).

ветствует форме правильных многогранников. По-

12.

L. Noirot, J. Nucl. Materials 447, 166 (2014).

казано, что коэффициенты диффузии для пузырь-

ков большинства рассмотренных размеров пропор-

13.

M. Veshchunov and V. Shestak, J. Nucl. Materials

циональны концентрации междоузельных атомов в

376, 174 (2008).

ОЦК-матрице.

14.

L. Verma, L. Noirot, and P. Maugis, J. Nucl.

Построена кинетическая модель, описывающая

Materials 528, 151874 (2020).

эволюцию функции распределения пузырьков по

размерам. Решение этой модели для времени в

15.

J. Evans, J. Nucl. Materials 210, 21 (1994).

1 мкс, что в 10 раз больше времени рассмотрен-

16.

А. С. Антропов, В. Д. Озрин, В. В. Стегайлов, В.

ных МД-расчетов, показывает длительное сохране-

И. Тарасов, ЖЭТФ 156, 125 (2019).

ние популяции малоподвижных пузырьков в матри-

17.

A. Antropov and V. Stegailov, J. Nucl. Materials 533,

це.

152110 (2020).

Расчеты выполнены с использованием Супер-

компьютера “Десмос” СКЦ ОИВТ РАН.

18.

A. Antropov and V. Stegailov, J. Nucl. Materials 551,

152942 (2021).

Финансирование. Исследование осуществлено

при частичной финансовой поддержке в рамках

19.

E. Gruber, J. Appl. Phys. 38, 243 (1967).

Программы фундаментальных исследований НИУ

20.

S. Chandrasekhar, Rev. Mod. Phys. 15, 59 (1943).

ВШЭ.

21.

E. Moore, L. R. Corrales, T. Desai, and R. Devana-

than, J. Nucl. Materials 419, 140 (2011).

ЛИТЕРАТУРА

22.

S. Murphy, A. Chartier, L. Van Brutzel, and J.-

P. Crocombette, Phys. Rev. B 85, 144102 (2012).

M. Tonks, D. Andersson, R. Devanathan, R. Dubo-

urg, A. El-Azab, M. Freyss, F. Iglesias, K. Kulacsy,

23.

A. Jelea, R.-M. Pellenq, and F. Ribeiro, J. Nucl.

G. Pastore, S. R. Phillpot et al., J. Nucl. Materials

Materials 444, 153 (2014).

504, 300 (2018).

24.

X.-Y. Liu and D. Andersson, J. Nucl. Materials 462,

J. Rest, M. Cooper, J. Spino, J. Turnbull, P. Van

8 2015

Uffelen, and C. Walker, J. Nucl. Materials 513, 310

(2019).

25.

L. Yang and B. Wirth, J. Nucl. Materials 544, 152730

(2021).

H. Trinkaus and B. Singh, J. Nucl. Materials 323,

26.

Z. Xiao, Y. Wang, S. Hu, Y. Li, and S.-Q. Shi, Comp.

229 (2003).

Mat. Sci. 184, 109867 (2020).

G. Pastore, D. Pizzocri, C. Rabiti, T. Barani,

27.

K. Nogita and K. Une, Nucl. Inst. and Met. in Phys.

P. Van Uffelen, and L. Luzzi, J. Nucl. Materials 509,

Res. Sec. B: Beam Interact. with Mat. and Atoms

687 (2018).

141, 481 1998

Z. Qian, W. Liu, R. Yu, Y. Tao, D. Yun, and L. Gu,

28.

G. Greenwood, A. Foreman, and D. Rimmer, J. Nucl.

J. Nucl. Materials 556, 53188 (2021).

Materials 1, 305 1959

D. Olander and D. Wongsawaeng, J. Nucl. Materials

29.

H. Xiao, C. Long, X. Tian, and S. Li, Mat. and Design

354, 94 (2006).

74, 55 2015

M. Veshchunov, V. Ozrin, V. Shestak, V. Tarasov,

30.

S. Hu, W. Setyawan, V. V. Joshi, and C. A. Lavender,

R. Dubourg, and G. Nicaise, Nucl. Eng. and Design

J. Nucl. Materials 490, 49 2017

236, 179 (2006).

31.

D. Yun, M. A. Kirk, P. M. Baldo, J. Rest, A. M. Ya-

A. Volkov and A. Ryazanov, J. Nucl. Materials 273,

cout, and Z. Z. Insepov, J. Nucl. Materials 437, 240

155 (1999).

2013

R. E. Voskoboinikov and A. E. Volkov, J. Nucl.

32.

D. Yun, J. Rest, W. Zhang, X. Xie, W. Liu, and

Materials 282, 66 (2000).

L. Gu, J. Nucl. Materials 540, 152409 (2020).

10.

R. E. Voskoboinikov and A. E. Volkov, J. Nucl.

33.

W. Zhang, D. Yun, and W. Liu, Materials 12, 2354

Materials 297, 262 (2001).

(2019).

212

ЖЭТФ, том 163, вып. 2, 2023

34.

G. Smirnov and V. Stegailov, J. Phys.: Cond. Mat.

45.

A. Shamsutdinov, M. Khalilov, T. Ismagilov, A. Pi-

31, 235704 (2019).

ryugin, S. Biryukov, V. Stegailov, and A. Timo-

feev, in International Conference on Mathematical

35.

D. E. Smirnova, A. Y. Kuksin, S. V. Starikov,

Modeling and Supercomputer Technologies (Springer,

V. V. Stegailov, Z. Insepov, J. Rest, and A. M. Ya-

2020), p. 401.

cout, Modell. Simul. Mater. Sci. Eng. 21, 035011 2013

46.

W. Jäger, R. Manzke, H. Trinkaus, G. Crecelius,

36.

D. Yun, Y. Miao, R. Xu, Z. Mei, K. Mo, W. Moha-

R. Zeller, J. Fink, and H. Bay, J. Nucl. Materials

med, B. Ye, M. J. Pellin, and A. M. Yacout, J. Nucl.

111, 674 (1982).

Materials 471, 272 2016

47.

M. Methfessel, D. Hennig, and M. Scheffler, Phys.

37.

B. Beeler, S. Hu, Y. Zhang, and Y. Gao, J. Nucl.

Rev. B 46, 4816 (1992).

Materials 530, 151961 (2020).

48.

W. Tyson, Canadian Metal. Quart. 14, 307 (1975).

38.

B. Beeler, M. W. Cooper, Z.-G. Mei, D. Schwen, and

Y. Zhang, J. Nucl. Materials 543, 152568 (2021).

49.

Z. Insepov, J. Rest, A. Yacout, A. Y. Kuksin, G. Nor-

man, V. Stegailov, S. Starikov, and A. Yanilkin,

39.

J. French and X.-M. Bai, J. Nucl. Materials 565,

J. Nucl. Materials 425, 41 (2012).

153744 (2022).

40.

I. Novoselov, A. Yanilkin, A. Shapeev, and E. Pod-

50.

V. V. Dremov, P. V. Chirkov, and A. V. Karavaev,

ryabinkin, Comp. Mat. Sci. 164, 46 (2019).

Sci. Rep. 11, 1 (2021).

41.

R. Ryltsev and N. Chtchelkatchev, J. Mol. Liq. 349,

51.

D. Seitov, K. Nekrasov, A. Y. Kupryazhkin, S. Gupta,

118181 (2022).

and A. Usseinov, Phys. Res. Sec. B 476, 26 (2020).

42.

A. Stukowski, Modell. Simul. Mater. Sci. Eng. 18,

52.

N. Kondratyuk, V. Nikolskiy, D. Pavlov, and V. Ste-

015012 (2010).

gailov, Int. J. High Perf. Comp. Appl. 35, 312 (2021).

43.

Г. Э. Норман, В. В. Стегайлов, Мат. Мод. 24, 3

53.

A. Leenaers, W. Van Renterghem, and S. Van den

(2012).

Berghe, J. Nucl. Materials 476, 218 (2016).

44.

A. P. Thompson, H. M. Aktulga, R. Berger,

54.

H. Xie, N. Gao, K. Xu, G.-H. Lu, T. Yu, and F. Yin,

D. S. Bolintineanu, W. M. Brown, P. S. Crozier,

Acta Materialia 141, 10 (2017).

P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Ngu-

yen et al., Computer Physics Communications 271,

55.

A. V. Nazarov, A. A. Mikheev, and A. P. Melnikov,

108171 (2022).

J. Nucl. Materials 532, 152067 (2020).

213