Квазиодномерная решеточная модель топливной ячейки Текст научной статьи по специальности «Физика»

38 ТРУДЫ БГТУ. 2015. № 6. Физико-математические науки и информатика. С. 38-42

УДК 531.19

Я. Г. Грода

Белорусский государственный технологический университет

КВАЗИОДНОМЕРНАЯ РЕШЕТОЧНАЯ МОДЕЛЬ ТОПЛИВНОЙ ЯЧЕЙКИ

Предложена одномерная решеточная модель топливной ячейки на основе иттрий-стабили-зированного диоксида циркония, допускающая численное моделирование процесса переноса заряда в рамках динамического метода Монте-Карло.

На основании разработанной решеточной модели топливной ячейки выполнено численное моделирование процесса переноса заряда в рамках динамического метода Монте-Карло.

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

Ключевые слова: топливная ячейка, решеточная модель, динамический метод Монте-Карло, вероятность перехода, энергия активации.

Ya. G. Groda

Belarusian State Technological University

QUASI-ONE-DIMENSIONAL LATTICE MODEL OF А FUEL CELL

The one dimensional lattice model of the fuel cell on the basis of Yttria-Stabilized Zirconia is proposed.

On the basis of the lattice model of the fuel cell in frame of kinetic Monte Carlo algorithm the numerical simulation of the charge transfer is realized.

It is shown that the electrostatic interaction between ions makes a significant contribution to the activation energy of migration of particles. On the other hand, the fluctuations of the energy barriers increase the activation energy only slightly. It is found that when electrodes are blocked electrical double layers are formed in the electrode areas. The thickness of the electrical double layer is not more than a few lattice constants.

Key words: fuel cell, lattice model, kinetic Monte Carlo, transition probability, activation energy.

Введение. В последнее время в литературе широко обсуждается моделирование твердотельных топливных элементов на макроскопическом уровне [1]. При этом основное внимание обращается на рассмотрение процесса переноса заряда через электроды элемента и изучение интегральных физико-химических процессов в самом элементе.

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

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

происходит на временах, значительно больших наносекундного диапазона. Например, на поверхности диоксида циркония, стабилизированного иттрием при 750 К, поверхностная диффузия иона кислорода на длине 8 А (примерно равной среднему расстоянию между дефектами решетки) происходит за 1 мк, а процесс поглощения одного молекулярного кислорода при давлении 0,01 атм - примерно за 0,5 мк. Кроме того, при МД-моделировании необходимо обрабатывать чрезвычайно большое число молекул окружающей среды, в которой находится изучаемая система. Поэтому наиболее выгодным с точки зрения временных затрат на моделирование является использование для таких систем динамического метода Монте-Карло (ДМК). При указанном подходе время измеряется в абстрактных шагах алгоритма Монте-Карло (МКШ). Это позволяет рассмотреть такие явления, как адсорбция или диффузия за один МКШ в отличие от МД-моделирования, в котором используются временные шаги малой величины, в результате чего отслеживаются промежуточные состояния атомов. Физическое время, соответствующее

КвазиоАномерная решеточная модель топливной ячейки

39

1 МКШ, связано с вероятностью, с которой происходят случайно выбранные события. Основное достоинство ДМК состоит в том, что данный подход позволяет сосредоточить внимание на сложностях, обусловленных взаимодействием множества отдельных процессов.

Квазиодномерная решеточная модель топливной ячейки. Одной из простейших моделей изучаемой системы, допускающей аналитическое рассмотрение, является модель Гуи - Чепмена. В рамках данного подхода предполагается, что поток ионов кислорода между катодом и анодом может быть описан соотношением

или

J = -Д ^ + ,

дх

_ дс оЕ J = -Д^ — +-,

Эх

Ч

(1)

(2)

где Дг и ВСь - коэффициенты кинетической и химической диффузии соответственно; В -абсолютная подвижность иона; ч - заряд иона; Е - электрическое поле; о - проводимость.

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

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

На рис. 1 представлена упрощенная одномерная дискретная модель активированного электролита, соответствующая модели Гуи -Чепмена, в которой весь массив электролита разбит на вертикальные слои, содержащие по три ячейки, которые могут быть заняты ионом или вакансией [2].

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

О о о о о о о о

• • • • • • • •

о о о о о о о о

Рис. 1. Одномерная решеточная модель топливного

элемента на основе иттрий-стабилизированного циркония. Темными точками показаны положительные ионы, серыми - отрицательные, белыми - вакансии.

При МК-моделировании подвижными являются только отрицательно заряженные ионы, совершающие переходы в соседние (слева либо справа) ячейки

Вычисление энергии активации в рамках динамического метода Монте-Карло. Предложенная решеточная модель является удобной для численного моделирования в рамках динамического метода Монте-Карло (ДМК). Основное отличие ДМК от классического алгоритма Метрополиса [3, 4] состоит в явном вычислении энергии активации еа, которая наряду с вероятностью перехода частицы на соседний слой является ключевым понятием данного подхода.

При объемной диффузии вероятность перехода частицы в соседнее положение определяется выражением аррениусовского вида

(

Р = Р0ехр

Л

квТ

(3)

где Р0 - нормировочный предэкспоненциаль-ный фактор; кв - постоянная Больцмана; Т -температура.

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

(4)

£ = £0 + 5е,

где е0 - средняя величина межузельного барьера в системе; 5е - некоторая добавка, которая является случайной величиной, подчиняющейся заданному распределению, например равномерному, гауссовому или экспоненциальному, либо зависит от положения ячейки и от времени.

Соотношение (4) позволяет рассматривать предлагаемую модель как один из вариантов изученной ранее неупорядоченной решетки [5].

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

поле, то связанная с этим движением работа равна:

£ Р = а0 q(U + Еюп + ЕаВ X (5)

где а0 - расстояние, на которое смещается ион; U - приложенная к электродам разность потенциалов; Еюп - вклад в общее электрическое поле со стороны подвижных и неподвижных ионов электролита; Есец - вклад в общее электрическое поле со стороны ионов, расположенных на том же слое, что и движущийся ион.

Фактически работа (5) приводит к изменению относительной высоты энергетического барьера, который необходимо преодолеть частице при переходе из одной ячейки в соседнюю, но в рассмотренном статическом приближении детали этого изменения не исследуются. Формально они могут быть учтены путем введения симетризационного множителя а. В дальнейшем мы будем полагать а = 0,5, что соответствует малому отличию в высоте энергетического барьера при прямом и обратном прыжке.

Таким образом, суммарная энергия активации, необходимая для перехода иона из ячейки j в соседнюю k, будет состоять из энергии активации объемной диффузии (4) и вклада электрического поля 2a(k - j)eP. Это позволяет записать вероятность перехода в следующем виде:

(

P = P0exp

£0 + б£ - 2a(k - j)£

kBf

Л

(6)

Величина электрического поля, обусловленного воздействием ионов электролита на слой у, вычисляется путем суммирования зарядов всех слоев, кроме слояу. При этом суммируются как положительно, так и отрицательно заряженные ионы. Суммарный положительный заряд слоя к вносит положительный вклад в поле Еоп, если к <у, и отрицательный, если к < у. Следовательно, для поля Еоп получаем:

Егюп (j) = -1 2£

N

I Q(k) - I Q(k)

k=j+i

k=i

(7)

где 8 - диэлектрическая проницаемость электролита; Q(к) - заряд к-го слоя.

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

Есе11 (j)=±-2£(Q< j)-q),

(8)

где положительный знак соответствует перемещению отрицательно заряженного иона на слой у + 1, а отрицательный - на слой у - 1.

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

В случае равномерного распределения высот межузельных барьеров на интервале [е0 - 5е; е0 + 5е] такой переход является тривиальным и определяется соотношением

ег =£0 + 25е(хг - 0,5). (9)

Параметры модели. Для моделирования движения ионов кислорода по решетке используются следующие физические параметры. Внешнее электрическое поле принято равным и = = 1 В/мкм = 106 В/м. Параметр решетки а0 = = 0,737 нм принят согласно данным работы [2]. В результате коэффициент взаимодействия двух-зарядных ионов кислорода на ближайших узлах решетки 2е2 / 80еа0 = 0,784 • 10-19 Дж, где е - заряд электрона, е0 - электрическая постоянная, а диэлектрическая постоянная среды е = 100 принята с учетом частичной экранировки ионов электронной подсистемой твердого электролита. После деления на тепловую энергию квТ при Т = 1000 К получим коэффициент при количестве частиц в уравнении (6) с учетом соотношения (8), равный 0,57 / t, где t = Т/ 1000. Влияние внешнего поля определяется коэффициентом 2еЕа0 / квТ = 0,017 / t.

Моделируемая решетка содержит 210 слоев. Сама процедура моделирования состоит из 1,01 • 106 МКШ, ' из которых первые 104 шагов отводятся для перехода системы в равновесное состояние (эквилибризации системы) и не учитываются при последующем усреднении результатов моделирования. При этом после выполнения каждых 100 шагов производится перестройка решетки, заключающаяся в новой расстановке случайных межузельных барьеров в соответствии с выбранным распределением.

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

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

41

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

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

Ситуация становится принципиально иной при высоких температурах. В этом случае можно утверждать, что из-за относительно высокой подвижности частиц заряд слоя меняется хаотически при переходе от слоя к слою. При этом его среднее значение оказывается близким к нулю.

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

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

На рис. 2 представлена зависимость от обратной температуры логарифма суммарного заряда в случае упорядоченной решетки с меж-узельным барьером е0 и при различных значениях параметра 5е, определяющего максимально возможное случайное отклонение величины межузельного барьера от среднего значения.

Необходимо отметить, что для оптимизации временных затрат и повышения точности результатов моделирование при различных параметрах 5е велось при разных значениях нормировочной константы Р0. Ее значения, как и отмечалось выше, определялись в ходе серии пробных численных экспериментов и принимались равными 1,445е0 для упорядоченной решетки, 1,87е0 для системы с 5е = = 0,5е0 и 2,55е0 для системы с 5е = 1,0е0. Представленные на рис. 2 данные приведены к единой системе отсчета, отвечающей упорядоченной решетке.

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

ределить среднюю энергию активации рассматриваемой модели.

0,8 1,2 1,6 2,0 2,4 2,8 Обратная температура, 1 ООО / Т

Рис. 2. Зависимость от обратной температуры логарифма суммарного заряда, прошедшего через электрод, в случае упорядоченной решетки

(кривая 1), при 5е = 0,5 (2) и 5е = 1,0 (3). Точками представлены результаты моделирования, сплошными линиями - результаты линейной аппроксимации данных моделирования

На рис. 3 представлена зависимость средней энергии активации от величины параметра 5е.

0,2 0,4 0,6 0,8 1,0

Коэффициент масштаба равномерного распределения, 5е / е0

Рис. 3. Зависимость средней энергии активации от величины коэффициента масштаба равномерного распределения

Анализируя рис. 3, можно отметить, что с увеличением параметра 5е, т. е. с ростом максимально возможного отклонения величины межузельного барьера от среднего значения, энергия активации увеличивается. При этом следует подчеркнуть, что указанный рост энергии активации происходит приблизительно по линейному закону вида

£a §£ шающей среднюю. В случае квазиодномерной

£" =2,752 +0,673—. (10) системы эти барьеры фактически «разрезают»

0 0 исходную решетку на отдельные части, и пере-

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

увеличении 5e возрастает число барьеров с ненным, что и проявляется через увеличение

энергетической высотой, существенно превы- средней энергии активации.

Литература

1. Zhu H., Kee R. J. A general mathematical model for analyzing the performance of fuel-cell membrane-electrode assemblies // J. Power Source. 2003. Vol. 117. P. 61-74.

2. 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.

3. Equation of state calculation by fast computing machines / N. Metropolis [et al.] // The Journal of Chemical Physics. 1953. Vol. 21, no. 6. P. 1087-1092.

4. Uebing C., Gomer R. A Monte Carlo study of surface diffusion coefficients in the presence of ad-sorbate-adsorbate interactions // The Journal of Chemical Physics. 1991. Vol. 95, no. 10. P. 7626-7652.

5. Вихренко В. С., Грода Я. Г., Бокун Г. С. Равновесные и диффузионные характеристики интеркаляционных систем на основе решеточных моделей. Минск: БГТУ, 2008. 326 с.

References

1. Zhu H., Kee R. J. A general mathematical model for analyzing the performance of fuel-cell membrane-electrode assemblies. J. Power Source, 2003, vol. 117, pp. 61-74.

2. 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.

3. Metropolis N., Rosenbluth A. W., Marshall N., Rosenbluth M. N., Teller A. H., Teller E. Equation of state calculation by fast computing machines. The Journal of Chemical Physics, 1953, vol. 21, no. 6, pp.1087-1092.

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. Vikhrenko V. S., Groda Ya. G., Bokun G. S. Ravnovesnyye i diffuzionnyye kharakteristiki interkal-yatsionnykh sistem na osnove reshetochnykh modeley [Equilibrium and diffusional properties of the intercalation systems in frame of lattice models]. Minsk, BSTU Publ., 2008. 326 p.

Информация об авторе

Грода Ярослав Геннадьевич - кандидат физико-математических наук, доцент, заведующий кафедрой теоретической механики. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: [email protected]

Information about the author

Groda Yaroslav Gennad'yevich - Ph. D. (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]

Поступила 12.03.2015