ВАКАНСИОННАЯ МИГРАЦИЯ КАТИОНОВ В ДИОКСИДЕ УРАНА. МОЛЕКУЛЯРНО-ДИНАМИЧЕСКОЕ МОДЕЛИРОВАНИЕ

Диоксид урана UO2 широко используется в качестве топлива для ядерной энергетики. Такие его свойства, как отсутствие фазовых переходов и взаимодействия с цирконием, ниобием, нержавеющей сталью при высоких температурах, позволяют эффективно применять его в ядерных реакторах, получая высокий КПД. Этими факторами и обусловлен выбор диоксида урана в качестве исследуемого модельного вещества. Диоксид урана, как и другие оксиды урана, используется также как промежуточный продукт при производстве других урановых соединений, главным образом фторидов. Все оксиды урана являются наиболее устойчивыми его соединениями, в связи с чем широко используются как для хранения урана, так и как промежуточное звено между урановорудным и фторидными урановыми производствами.

Точечные дефекты в кристаллических решетках твердых веществ играют важную роль в миграции атомов. Они влияют на механические свойства материала, его проводимость, подвижность носителей заряда и др. Образование и миграция дефектов типа «вакансия» напрямую связано с явлением самодиффузии веществ

Создание модели физических объектов в большей степени помогает избавиться от лишних временных, финансовых или кадровых затрат на изучение самого объекта. Понимание основных характеристик и функций объекта, цель исследования, а также корректная расстановка приоритетов между исследуемыми процессами деятельности объекта помогает создать оптимальную по всем пунктам модель. В нашей работе такой моделью является компьютерная программа, описывающая поведение множества материальных точек, связанных между собой физическими законами. Объектом, по которой построена модель, является трехмерная кристаллическая решетка диоксида урана. Удалив из модели одну или более специальным образом выбранные точки, можно получить аналог такого физического явления, как образование точечного дефекта в кристалле твердого тела вакансионного типа. Необходимая производительность вычислений достигнута с использованием программно-аппаратной платформы CUDA компании NVIDIA, обеспечившей возможность организации доступа к набору инструкций графического ускорителя и управления его памятью при организации параллельных вычислений

В основе метода молекулярной динамики (МД) лежит модельное представление о многоатомной молекулярной системе, в которой все атомы представлены материальными точками, движение которых описывается в классическом случае уравнениями Ньютона

1) имеет массу, радиус-вектор и скорость соответственно , , , где ,

2) взаимодействует с остальными посредством сил , где потенциальная энергия взаимодействия системы из N частиц,

3) взаимодействует с внешними полями посредством силы .

Тогда эволюция данной модели будет описываться системой 2N обыкновенных дифференциальных уравнений движения.

В работе задействуется один из наиболее используемых потенциалов в методе молекулярной динамики – потенциал Букингема:

(1)

Потенциал (1) описывает отталкивание электронных оболочек в экспоненциальной форме, которая является наиболее теоретически обоснованной из простых функциональных зависимостей. Он не учитывает затухания дисперсионного притяжения на малых расстояниях порядка суммы радиусов взаимодействующих частиц, на которых не справедливо мультипольное разложение

Заключается алгоритм МД в том, что кристалл представляется системой недеформируемых, положительных и отрицательных, ионов, эволюционирующей во времени. Ионы перемещаются по законам Ньютона, а силы взаимодействия определяются парными потенциалами Uij(Rij)

Потенциалы Uij(Rij) имеют общий вид

(2)

где QiQj – заряды i-го и j-го ионов;

Rij – расстояние между этими ионами;

KE – константа закона Кулона;

UijОб (Rij) – некулоновский потенциал взаимодействия электронных оболочек (1).

Потенциалы UijОб (Rij) представляют в различных формах, обычно включающих в себя слагаемые, моделирующие отталкивание перекрывающихся электронных оболочек на малых расстояниях и дисперсионное притяжение на сравнительно больших. Наиболее известными из этих форм являются потенциал Букингема (1) и потенциал Леннарда-Джонса.

В формуле (1) А и В – константы, характеризующие отталкивание оболочек, C6, C8 – константы, описывающие дипольную и квадрупольную составляющие дисперсионного притяжения.

Перед началом расчёта ионам присваиваются некоторые начальные координаты и скорости (например – координаты, соответствующие узлам идеальной кристаллической решётки моделируемого соединения и скорости, соответствующие Максвелловскому распределению при заданной температуре), после чего начинается численное пошаговое интегрирование уравнений движения ионов во времени. На каждом k-ом шаге производятся следующие действия:

1) рассчитываются действующие на каждый ион силы

(3)

2) вычисляются новые скорости и новые координаты ионов

(4)

(5)

Формулы (3-5) справедливы при нулевых граничных условиях (конечный кристаллит из N частиц в вакууме) без компенсации перемещения, вращения и дрейфа температуры, возникающих из-за вычислительных погрешностей (алгоритмы компенсации приведены ниже). Однако этих формул достаточно, чтобы показать возможность эффективного распараллеливания по схеме SIMD: очевидно, что основные этапы алгоритма заключаются в проведении над каждым из ионов поочерёдно одних и тех же операций.

Наиболее критичным участком алгоритма является расчёт результирующих сил, действующих на каждый из ионов со стороны остальных. Этот расчёт необходим на каждом шаге молекулярной динамики, а его объём квадратичен по количеству частиц N, так как для каждой из N частиц необходимо выполнить суммирование сил, действующих со стороны N-1 остальных:

(6)

Поскольку для реалистичного моделирования модельные кристаллиты должны содержать десятки и сотни тысяч ионов, объём расчёта (6) очень велик по сравнению с расчётами на остальных этапах алгоритма. Именно этот расчёт имеет смысл в первую очередь реализовать на графических процессорах. Нам такая реализация позволила на порядки ускорить молекулярно-динамическое моделирование ионных кристаллов (диоксида урана)

В какой-то момент атом может получить от соседей такой избыток энергии, что он займёт соседнее положение в решётке. Так осуществляется миграция (перемещение) точечных дефектов в объёме кристаллов

Если один из атомов, окружающих вакансию, переместится в вакантный узел, то вакансия соответственно переместится на его место. Последовательные элементарные акты перемещения определённой вакансии осуществляются разными атомами. На Рисунке 1 показано, что в слое плотноупакованных шаров (атомов) для перемещения одного из шаров в вакантное место он должен раздвинуть соседние шары.

Рисунок 1 - Точечные дефекты

Примечание: по ист. [7]

DOI:10.60797/IRJ.2024.143.113.1

Следовательно, для перехода из положения в узле, где энергия атома минимальна, в соседний вакантный узел, где энергия также минимальна, атом должен пройти через состояние с повышенной потенциальной энергией, преодолеть энергетический барьер. Для этого и необходимо атому получить от соседей избыток энергии, который он теряет, «протискиваясь» в новое положение

[8]

.

Моделирование миграции катионов в настоящей работе включало следующие этапы:

1) создание модель кристаллической решетки диоксида урана;

2) удаление одной молекулы UO2 из решетки диоксида урана для формирования вакансии;

3) моделирование поведение кристаллической решетки в течение некоторого интервала времени,

Кристаллическая решетка диоксида урана имеет гранецентрированную кубическую форму. Теоретическая плотность – 10,96 г/см3, температура плавления 2875±45 °С

В качестве исследуемой модели используется нанокристалл кубической формы размером 22×22×22 Å (ангстрем), содержащий 256 атомов 92238U и 512 атомов 816O (всего 768 атомов) (рисунок 2). Потенциалом взаимодействия выбран потенциал Букингема (1), дополненный слагаемым, характеризующим кулоновское отталкивание частиц (2).

В качестве граничных условий используется математический аналог изолирования кристаллита в вакууме. Кроме того, используется термостат Берендсена, представляющий собой алгоритм для изменения масштаба скоростей частиц в моделировании молекулярной динамики для управления температурой моделирования. В связи с конечностью кристаллической решетки, атомы диоксида урана начинают приобретать нефизичные значения скоростей. Термостат позволяет поддерживать скорости частиц в определенных пределах, отводя лишнюю энергию от системы.

Рисунок 2 - Кристаллическая решетка диоксида урана в начальный момент времени: 

серые большие шары – атомы урана; красные меньшие – атомы кислорода

DOI:10.60797/IRJ.2024.143.113.2

После создания модели решетки в кристалл вносили вакансию, состоящую из двух атомов кислорода и атома урана, т. е. удалить из кристалла одну из молекул (рисунок 3). Вакансию размещали в центре кубической решетки, чтобы при диффузии частиц она как можно дольше находилась внутри нанокристалла.

После релаксации кристаллической решетки в течение достаточного количества времени исследовали поведение катионов их миграцию внутри кристалла. При этом, удостоверялись в том, что вакансия находилась внутри исследуемой зоны кристаллической решетки. Для этого был создан метод регистрации вакансии в любой момент времени, основанный на понятии координационной сферы. Координационная сфера – совокупность атомов или молекул, расположенных на одинаковом расстоянии от рассматриваемой частицы (Рисунок 4). Метод регистрации катионной вакансии состоит в следующем.

Рисунок 3 - Наглядное изображение вакансии внутри кристаллической решетки

DOI:10.60797/IRJ.2024.143.113.3

Рисунок 4 - Пример координационной сферы

DOI:10.60797/IRJ.2024.143.113.4

В построенной кристаллической решетке выделяется некоторая область, внутри которой находится предполагаемая вакансия. Поочередно выбираются катионы (атомы урана) из этой области и определяется их первая координационная сфера, т. е. находятся «ближайшие соседи». Из-за того, что одного урана не хватает, у некоторых катионов вместо положенного числа «соседей» будет на одного меньше. Тогда очевидно, что эти катионы будут расположены вблизи самой вакансии (вокруг нее). Остается лишь из известных координат катионов найти усредненные, как раз и являющиеся координатами вакансии (рисунок 5).

Рисунок 5 - Метод регистрации катионной вакансии

DOI:10.60797/IRJ.2024.143.113.5

В первую очередь необходимо установить наличие вакансии в кристалле вышеприведенным методом и оценить факт ее перемещения внутри кристалла с течением времени. По истечении времени моделирования были получены координаты как самих атомов урана и кислорода, так и координаты вакансии в каждый момент времени. На Рисунке 6 приведен вид кристаллической решетки в последний момент времени расчета, где фиолетовой сферой обозначено местоположение катионной вакансии.

Таким образом, факт нахождения самой вакансии в кристалле подтвержден. Далее следует установить наличие миграций данной вакансии и катионов внутри кристалла с течением времени. Для этого вычисляли квадрат смещения σ2 вакансии и атомов вокруг нее относительно их начального положения в каждый момент времени:

(8)

где X0, Y0, Z0 – начальные координаты вакансии и катионов;

Xi, Yi, Zi – координаты вакансии и катионов в текущий момент времени.

Получив нужные значения, строили график зависимости квадрата смещения вакансии от номера итерации для двух значений температур (Рисунок 7).

Из полученного графика зависимости видно, что вакансия не остается в первоначальном положении, но и не перемещается на расстояние большее, чем постоянная решетки (5,47 Å) и периодически возвращается в исходное местоположение. В связи с этим, факт наличия миграции вакансии вокруг своего первоначального положения считается подтвержденным. Кроме того, очевидно, что при температуре кристалла 3000 К вакансия смещается на большее расстояние от первоначального положения.

Рисунок 6 - Кристаллическая решетка диоксида урана в конечный момент времени моделирования:

большие шары – атомы урана; меньшие – атомы кислорода; фиолетовые – местоположение вакансии

DOI:10.60797/IRJ.2024.143.113.6

Рисунок 7 - Зависимость квадрата смещения вакансии от времени при температурах 2000 К и 2800 К

DOI:10.60797/IRJ.2024.143.113.7

Проследить за траекторией движения вакансии может помочь построение двумерных графиков проекции координаты вакансии на координатные плоскости xOy, zOy и zOx (Рисунки 8-10), из которых так же видно, что вакансия не перемещалась от своего первоначального положения дальше 3 Å.

Рисунок 8 - Проекция координат катионной вакансии на плоскость xOy

DOI:10.60797/IRJ.2024.143.113.8

Рисунок 9 - Проекция координат катионной вакансии на плоскость zOx

DOI:10.60797/IRJ.2024.143.113.9

Из проекций видно, что при температуре в 3000 К вакансия может перемещаться на большие расстояния, чем при 2000 К. Характер миграции вакансии носит сугубо вероятностный характер, т. е. повторное моделирование при тех же самых условиях может дать абсолютно другие изображения проекций.

Иное дело обстоит с миграцией катионов вблизи вакансии. На Рисунке 11 приведена зависимость среднего квадрата смещения атомов урана от номера итерации для двух значений температур. В отличие от вакансии катионы не колеблются вблизи своего исходного положения, и со временем расчета повышается их величина смещения, что является фактом наличия миграции катионов.

Рисунок 10 - Рост среднего квадрата смещения катионов со временем расчета

DOI:10.60797/IRJ.2024.143.113.10

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

(9)

где 〈a2〉 – средний квадрат смещения, полученный из рисунка 11;

t – время процесса диффузии.

Температурная зависимость коэффициента диффузии имеет вид:

(10)

где E – энергия эффективной активации диффузии,

k – постоянная Больцмана, k=1,38∙10-23 Дж/K.

После преобразований формулы (10) можно получить зависимость логарифма коэффициента диффузии от обратной температуры. Угловой коэффициент полученной прямой будет равняться значению энергии эффективной активации диффузии.

(11)

Для построения необходимой зависимости были рассчитаны значения коэффициента диффузии для некоторых значений температур в диапазоне 2200-2900 К. В таблице 1 приведены рассчитанные значения коэффициента диффузии, обратной температуры, а также промежуточные величины. На рисунке 12 изображена вышеупомянутая зависимость.

Рисунок 11 - Зависимость логарифма коэффициента диффузии от обратной температуры

Примечание: пунктирной линией изображена аппроксимирующая функция

DOI:10.60797/IRJ.2024.143.113.12

Сравнивая уравнение (11) с аппроксимирующей функцией, приведенной на рисунке 12, делается вывод, что , тогда .

Полученное значение энергии эффективной активации диффузии можно сравнить с экспериментальными и расчетными значениями для одиночной катионной вакансии – 1,5-3,0 эВ

Таким образом, в работе предложена модель, описывающая поведение кристаллической структуры с точечным дефектом, построенная на методе молекулярной динамики. Модельными системами были кубические гранецентрированные кристаллиты размером в 765 частиц с точечным дефектом вакансионного типа. Граничные условия соблюдаются благодаря использованию термостата Берендсена, который поддерживает нужную температуру, отводя лишнюю энергию от решетки.

В качестве оценки миграции катионов построена зависимость логарифма коэффициента диффузии катионов от обратной температуры кристаллической решетки, предварительно были найдены сами значения коэффициента диффузии с помощью известной связи со средним квадратом смещения атомов в кристалле. Абсолютное значение углового коэффициента полученной с помощью линейной аппроксимации прямой дает значение эффективной энергии активации диффузии. Данная величина, равная 1,10 эВ, оказывается близкой к экспериментальным и расчетным значениям других авторов для одиночной катионной вакансии, находящимся в диапазоне 1,5-3,0 эВ. Несколько заниженный результат можно объяснить наличием помимо катионной вакансии двух анионных вакансий. Весь вакансионный кластер обладает меньшей подвижностью по сравнению с одиночной катионной вакансией.