CC BY
12
46
Труды БГТУ, 2016, № 6, с. 46-50
УДК 531.19
Я. Г. Грода, В. С. Вихренко
Белорусский государственный технологический университет
ВЛИЯНИЕ ПАРАМЕТРОВ СЛУЧАЙНОГО РАСПРЕДЕЛЕНИЯ ВЫСОТ МЕЖУЗЕЛЬНЫХ БАРЬЕРОВ НА ЭНЕРГИЮ АКТИВАЦИИ ПРОВОДИМОСТИ В КВАЗИОДНОМЕРНОЙ РЕШЕТОЧНОЙ МОДЕЛИ ТОПЛИВНОЙ ЯЧЕЙКИ
Рассматривается квазиодномерная решеточная модель топливной ячейки на основе иттрий-стабилизированного диоксида циркония, в которой ионы кислорода совершают термоактивированные переходы в вакантные решеточные узлы, а ионы иттрия и циркония принимаются неподвижными. Активационный барьер, преодолеваемый частицей при переходе, полагается состоящим из случайного межузельного барьера, величина которого подчиняется распределению Гаусса с заданными средним значением и дисперсией, и электростатической части, обусловленной взаимодействием между ионами. Для численного моделирования переноса заряда в системе использовался динамический метод Монте-Карло. Показано, что электропроводность обратно пропорциональна температуре, что позволяет определить ее энергию активации. Электростатическое взаимодействие приводит к возрастанию энергии активации электропроводности, в то время как увеличение среднеквадратичного отклонения высот межузельных барьеров от среднего значения вызывает ее снижение. Установлено, что при блокированных электродах в приэлектродных областях образуются двойные электрические слои толщиной не более нескольких постоянных решетки, которая слабо зависит от температуры и дисперсии гауссовского распределения барьеров.
Ключевые слова: топливная ячейка, решеточная модель, динамический метод Монте-Карло, вероятность перехода, энергия активации.
Ya. G. Groda, V. S. Vikhrenko
Belarusian State Technological University
INFLUENCE OF RANDOM DISTRIBUTION PARAMETERS OF INTERSITE BARRIER HEIGHTS ON CONDUCTIVITY ACTIVATION ENERGY IN QUASI-ONE-DIMENSIONAL LATTICE MODEL OF FUEL CELL
The quasi-one-dimensional lattice model of an oxide fuel cell is considered. Oxygen ions can perform thermally activated jumps over vacant lattice sites while positive yttrium and zirconium ions are considered as immobile. The activation barrier for particle jumps consists of a host system contribution that can fluctuate with Gaussian distribution around some energy level and electrostatic interactions between ions. Kinetic Monte Carlo algorithm is used for numerical simulation of charge transfer in the system. It is shown that the electrical conductivity is proportional to inverse temperature that allows us to evaluate the conductivity activation energy. The electrostatic interactions increase the latter while the host barrier fluctuations decrease it. It is found that for blocked electrodes an electric double layer is formed near the electrodes. The thickness of the double layer is of a few lattice constants and weakly depends on the temperature and Gaussian distribution dispersion.
Key words: fuel cell, lattice model, kinetic Monte Carlo, transition probability, activation energy.
Введение. В работе [1] было отмечено, что применение динамического метода Монте-Карло (ДМК) при моделировании твердотельных топливных элементов оказывается значительно более выгодным с точки зрения временных затрат по сравнению с молекулярно-динамическим (МД) моделированием. Это обусловлено в первую очередь тем, что в рамках ДМК время измеряется в абстрактных шагах алгоритма Монте-Карло (МКШ). Это позволяет рассмотреть такие явления, как адсорбция или диффузия, за один МКШ в отличие от МД-моделирования, в котором используются временные шаги малой величины и отслеживаются промежуточные состояния атомов.
С целью применения ДМК к изучению электрофизических свойств топливной ячейки на основе иттрий-стабилизированного диоксида циркония в работе [1] была предложена квазиодномерная решеточная модель системы, представленная на рис. 1.
В рамках этой модели весь массив электролита разбит на укрупненные ячейки, каждая из которых содержит по одной вакансии в подре-шетке кислорода и по два иона иттрия, что обеспечивает электронейтральность системы. Ионы иттрия и циркония малоподвижны, тогда как ионы кислорода могут перемещаться в вакантные узлы. На рис. 1 укрупненные ячейки показаны вертикальными слоями, в каждом из кото-
Я Г. Грода, В. С. Вихренко
47
рых положительные ионы представлены одним темным кружком. В используемой ниже модификации этой модели предполагается, что в каждой ячейке имеется два узла, доступных ионам кислорода (серые кружки), что обеспечивает дополнительные пути их миграции и возможность перестановки ионов вдоль горизонтальной оси, отсутствующей в чисто одномерных моделях. В каждом слое могут находиться либо два иона кислорода, либо две вакансии (светлые кружки), либо одна вакансия и один ион кислорода. Ионы кислорода могут совершать термоактивированные, с учетом действия электрического поля, переходы из одного слоя в другой.
О о о о о о о о
• • • • • • • •
о о о о о о о о
Рис. 1. Квазиодномерная решеточная модель топливного элемента на основе иттрий-стабилизированного циркония. Темными кружками показаны положительные ионы, серыми - отрицательные, белыми - вакансии.
Отрицательно заряженные ионы могут переходить в соседние (слева либо справа) ячейки
Проведенное в работах [2, 3] моделирование описанных систем с равномерным и экспоненциальным распределением высот межузельных барьеров показало, что электростатическое взаимодействие между ионами вносит значительный вклад в энергию активации миграции частиц. В свою очередь, флуктуации энергетических барьеров влияют на энергию активации незначительно. Также было установлено, что при блокированных электродах в приэлектродных областях образуются двойные электрические слои толщиной не более нескольких постоянных решетки. В настоящей работе представлены результаты моделирования системы с нормальным (гауссовским) распределением высот межузель-ных барьеров и исследовано влияние дисперсии распределения на энергию активации системы.
Вычисление энергии активации в рамках динамического метода Монте-Карло. Как было отмечено ранее [2], основное отличие ДМК от классического алгоритма Метрополиса [4] состоит в явном вычислении энергии активации Еа, которая определяет вероятность перехода частицы в соседний слой и является ключевым понятием данного подхода.
При объемной диффузии вероятность перехода частицы в соседнее положение Р определяется выражением аррениусовского вида
Р = Р0-1 ехр
КГ
(1)
где Ро - нормировочный предэкспоненциаль-ный фактор; кв - постоянная Больцмана; Т -температура.
Поскольку мигрирующие частицы не являются электронейтральными, необходимо учитывать влияние на энергию активации как неоднородности кристаллического поля, вызванное, например, наличием межзеренных границ в электролите, так и электрических полей. Это позволяет представить энергию активации в следующем виде:
Еа = Ев (Ео, а) + ЕР
(2)
где величина межузельного барьера Ев представляет собой случайную величину, подчиняющуюся нормальному распределению с математическим ожиданием Ео и среднеквадратичным отклонением (дисперсией) о.
Вклад ЕР электрических полей в энергию активации определяется электростатическим взаимодействием ионов и может быть записан как:
Ер = к(и + Ейоп) + Е
(3)
где к = ±1 и определяется направлением прыжка отрицательно заряженного иона со слоя . на слой . - 1 либо . + 1 соответственно; и - энергия внешнего электрического поля; Е.юп - вклад в общее электрическое поле со стороны подвижных и неподвижных ионов электролита всех его слоев, за исключением слоя у; Е^Ьее1 -вклад в общее электрическое поле со стороны ионов, расположенных в слое т. е. в том же слое, что и движущийся ион [1, 3]:
Е.. =
.юп
ЛЕ
( N
1
Л
I в(к) - X в(к)
^ к=. +1 к =1
ЛЕ
=—((к) - д),
(4)
(5)
где ЛЕ - энергия взаимодействия двухзарядных ионов кислорода, находящихся в соседних плоскостях в расчете на один ион; Q(k) - заряд к-го слоя модели.
Алгоритм моделирования. Для моделирования рассматриваемой системы по методу Монте-Карло используется алгоритм, в котором случайным образом выбирается произвольный отрицательно заряженный ион, занимающий слой у. После этого определяется направление его возможного перехода. Если ячейка-приемник к оказывается занятой другим отрицательным ионом, то переход считается
невозможным, так как двойное занятие ячеек не допускается. Тем не менее попытка этого перехода учитывается. Если ячейка-приемник свободна, то переход считается происходящим с вероятностью, определяемой соотношением (1), в котором нормировочная вероятность Р0 подбирается из условия недопущения Р > 1 непосредственно в ходе моделирования. Для принятия либо отклонения рассматриваемого перехода генерируется случайное число Рг из диапазона [0; 1], которое сопоставляется с Р. Если Р > Рг, то переход считается произошедшим, если Р < Рг, то переход отклоняется.
Повторение описанной процедуры N раз, где N - число рассматриваемых слоев, формирует один шаг алгоритма Монте-Карло (МКШ).
В процессе моделирования некоторое начальное количество МКШ отводится для экви-либризации системы, т. е. перехода системы в равновесное состояние, и не учитывается при последующем анализе результатов моделирования. Для процедуры усреднения результатов, полученных по различным реализациям предложенной модельной системы, использовался алгоритм, в котором после некоторого числа МКШ происходит перестройка моделируемой системы, заключающаяся в задании нового случайного распределения значений энергий межузельных барьеров.
Для определения величины заряда, проходящего через анод, а следовательно, и тока, протекающего через электролит, применяются периодические граничные условия.
Если же отказаться от периодических граничных условий и запретить ионам покидать моделируемую систему, то путем усреднения зарядов в слоях по времени может быть определено равновесное пространственное распределение отрицательно заряженных ионов во внешнем поле.
При установлении величин межузельных барьеров Ев, входящих в соотношение (2), следует отметить, что с помощью встроенных стандартных средств языков программирования можно, как правило, выбрать лишь равномерно распределенное псевдослучайное число хг из диапазона [0; 1]. Следовательно, необходимо осуществить переход от равномерно распределенной случайной величины к случайным числам Ев, распределенным согласно заданной функции распределения.
В случае нормального распределения такой переход может быть выполнен с помощью преобразования Бокса - Мюллера [5]:
где ф, r - независимые случайные величины, равномерно распределенные на интервале (0; 1].
Это преобразование является точным, однако его использование предъявляет повышенные требования к качеству генератора псевдослучайных чисел. Для одномерной последовательности случайных нормально распределенных чисел стандартный генератор языка Fortran 90 оказывается вполне применимым. Дополнительно псевдослучайные числа могут быть сделаны «более случайными» в процессе моделирования путем периодического сброса генератора и введения нового затравочного параметра, основанного на текущем системном времени. В результате будут получены два независимых числа z1 и z2, удовлетворяющих стандартному нормальному распределению с математическим ожиданием 0 и дисперсией 1.
Случайная величина барьера, распределенная нормально с заданным математическим ожиданием E0 и среднеквадратичным отклонением о, находится по формуле
EB E0 '
Oz,
(7)
z1 = cos(2^)V -2ln r, z2 = sin(2^)V -2ln r,
(6)
где в качестве г выступают последовательно г и г2 для перехода в данный узел из двух разных узлов соседнего слоя, выбранных с одинаковыми вероятностями, и для каждого узла последовательно разыгрываются новые числа и г2.
Параметры модели. Параметр укрупненной решетки а0 и средняя величина межузель-ного барьера Е0 могут быть приняты, согласно данным работы [1], равными 0,737 нм и 0,83 эВ соответственно.
Напряженность внешнего электрического поля Е полагается равной 106 В/м, что приводит к следующему значению энергии и
и = 2еЕа0 = 2,36 -10-22 Дж -1,5 -10-3 эВ. (8)
Энергия взаимодействия двухзарядных ионов кислорода ДЕ в соседних слоях определяется в соответствии с соотношением
2 2
АЕ =-= 7,87 10-21 Дж - 5 10-2эВ, (9)
££0а0
где е - заряд электрона; £0 - электрическая постоянная, а диэлектрическая постоянная среды £ = 1000 принята с учетом частичной экранировки ионов электронной подсистемой твердого электролита.
Моделируемая решетка содержит 128 слоев. Сама процедура моделирования состоит из 107 МКШ. При этом дополнительные первые 2 • 104 шагов отводятся для перехода системы в равновесное состояние (эквилибризация сис-
Я. Г. Грода, В. С. Вихренко
49
темы) и не учитываются при последующем усреднении результатов моделирования. После выполнения каждых 1оо шагов производится перестройка решетки, заключающаяся в новой расстановке случайных межузельных барьеров в соответствии с выбранным распределением.
Результаты моделирования и их обсуждение. Рассмотрение «закрытой» системы, в которой невозможно движение мигрирующих частиц через границы, позволяет изучить равновесное распределение зарядов в электролите (рис. 2).
0,15
а
« 0,10
н о
« « 0,05
I 0,00
(и Л -0,05
О
-0,10
40 60 80 100 Номер слоя
Рис. 2. Распределение заряда в модели при о = 0,12 эВ и Т = 1000 К
Анализ данного распределения показывает, что вблизи каждого из электродов образуется двойной электрический слой. Также можно отметить, что суммарный заряд такого двойного слоя и его ширина практически не зависят ни от температуры, ни от флуктуаций высот межузельных барьеров.
При пользовании в алгоритме моделирования периодических граничных условий появляется возможность определить суммарный заряд, проходящий через электрод за все время моделирования, т. е. фактически ток, протекающий через электролит.
На рис. 3 представлена зависимость от обратной температуры логарифма заряда, прошедшего через электроды, при различных значениях среднеквадратичного отклонения. Полученные зависимости являются практически линейными, так что углы наклона полученных прямых позволяют установить энергию активации электропроводности рассматриваемой модели.
Необходимо отметить, что для оптимизации временных затрат и повышения точности результатов моделирование при различных значениях параметра о велось при разных значениях нормировочной константы Р0 или, что эквивалентно, при различных значениях средней величины межузельного барьера.
0.6 0,8 1,0 1,2 Обратная температура, 1000 /Т
Рис. 3. Зависимость логарифма суммарного заряда, прошедшего через электрод, от обратной температуры для системы с о = 0,01 эВ и Е0 = 0,4 эВ (линия 1), о = 0,06 эВ и Е0 = 0,5 эВ (линия 2), о = 0,09 эВ
и Е0 = 0,6 эВ (линия 3), о = 0,11 эВ и Е0 = 0,7 эВ (линия 4). Точками обозначены результаты моделирования, сплошными линиями -результаты линейной аппроксимации данных моделирования
На рис. 4 показана зависимость средней энергии активации от дисперсии о. Представленные здесь данные приведены к единой системе, соответствующей средней величине межузельного барьера, равной 0,83 эВ. Это значение на рис. 4 показано пунктирной линией.
о !).;>Л П.П.' 1.1.1-
Среднеквадратичное отклонение ст, эВ
Рис. 4. Зависимость средней энергии активации от дисперсии
Заключение. Электростатическое взаимодействие приводит к незначительному (около 0,05 эВ при о = 0) повышению энергии активации. Увеличение же дисперсии о, напротив, вызывает снижение энергетического барьера. При о ~ 0,065 эВ средняя энергия активации становится равной средней величине межузель-ного барьера.
Полученная зависимость энергии активации от параметра о может быть объяснена тем, что при его увеличении возрастает число межузельных
барьеров с высотой, значительно меньшей средней. Тем самым у частиц появляется возможность мигрировать по системе, преодолевая
преимущественно относительно низкие энергетические барьеры, что и проявляется в снижении средней энергии активации.
Литература
1. Modak A. U., Lusk M. T. Kinetic Monte Carlo simulation of a solid-oxide fuel cell: I. Open-circuit voltage and double layer structure // Solid State Ionics. 2005. Vol. 176. P. 2181-2191.
2. Грода Я. Г. Квазиодномерная решеточная модель топливной ячейки // Труды БГТУ. 2015. № 6: Физ.-мат. науки и информатика. С. 38-42.
3. Charge Distribution Around Nanoscale Nonhomogeneities in Solid State Ionics / G. S. Bokun [et al.] // Nanomaterials: Applications and Properties. 2015. Vol. 4, no. 1. Art. # 01PCSI07 (4 p.).
4. Uebing C., Gomer R. A Monte Carlo study of surface diffusion coefficients in the presence of adsorbate-adsorbate interactions // The Journal of Chemical Physics. 1991. Vol. 95, no. 10. P. 7626-7652.
5. Box G. E. P., Muller M. E. A note on the generation of random normal deviates // The Annals of Mathematical Statistics. 1958. Vol. 29. P. 610-611.
References
1. Modak A. U., Lusk M. T. Kinetic Monte Carlo simulation of a solid-oxide fuel cell: I. Open-circuit voltage and double layer structure. Solid State Ionics, 2005, vol. 176, pp. 2181-2191.
2. Groda Ya. G. Quasi-one-dimensional lattice model of а fuel cell. Trudy BGTU [Proceedings of BSTU], 2015, no. 6: Physical-mathematical sciences and informatics, pp. 38-42 (In Russian).
3. Bokun G. S., Groda Ya. G., Lasovsky R. N., Vikhrenko V. S. Charge distribution around nanoscale nonhomogeneities in solid state ionics. Nanomaterials: Applications and Properties, 2015, vol. 4, no. 1, ай. # 01PCSI07 (4 p.).
4. Uebing C., Gomer R. A Monte Carlo study of surface diffusion coefficients in the presence of adsorbate-adsorbate interactions. The Journal of Chemical Physics, 1991, vol. 95, no. 10, pp. 7626-7652.
5. Box G. E. P., Muller M. E. A note on the generation of random normal deviates. The Annals of Mathematical Statistics, 1958, vol. 29, pp. 610-611.
Информация об авторах
Грода Ярослав Геннадьевич - кандидат физико-математических наук, доцент, заведующий кафедрой теоретической механики. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: [email protected]
Вихренко Вячеслав Степанович — доктор физико-математических наук, профессор, профессор кафедры теоретической механики. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: [email protected]
Information about the authors
Groda Yaroslav Gennad'yevich — PhD (Physics and Mathematics), Assistant Professor, Head of the Department of Theoretical Mechanics. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: [email protected]
Vikhrenko Vyacheslav Stepanovich — DSc (Physics and Mathematics), Professor, Professor, the Department of Theoretical Mechanics. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: [email protected]
Поступила 10.03.2016
