CC BY
15
ТЕОРЕТИЧЕСКАЯ МЕХАНИКА
УДК 531.19
Р. Н. Ласовский, Д. В. Гапанюк, Я. Г. Грода
Белорусский государственный технологический университет
МОДЕЛИРОВАНИЕ ТРЕХМЕРНОГО ТВЕРДОТЕЛЬНОГО ЭЛЕКТРОЛИТА
СО СЛЯБ ГЕОМЕТРИЕЙ
Рассматривается трехмерная модель твердотельного электролита со сляб геометрией, содержащей зерно и межзеренную прослойку, описываемую слоями, которые характеризуются дополнительными межузловыми энергетическими барьерами. Взаимодействие между подвижными ионами принимается состоящим из дальнодействующего электростатического отталкивания и короткодействующего притяжения. Система находится между двумя противоположно заряженными электродами. Выполняется моделирование описанной системы по методу Монте-Карло. При этом кулоновская энергия определяется суммированием по методу Эвальда с добавлением дополнительных «вакуумных областей» с обеих сторон сляба. Исследуется профиль концентрации подвижных ионов при различных условиях. Сопоставление результатов с аналитическими вычислениями показывает их количественное расхождение в пределах 15%.
Ключевые слова: твердотельный электролит, межзеренная граница, электрод, метод Монте-Карло, суммирование Эвальда, электропроводность, профиль концентрации.
R. N. Lasovsky, D. V. Gapanjuk, Ya. G. Groda
Belarusian State Technological University
MODELING OF THREE-DIMENSIONAL SOLID ELECTROLYTE WITH A SLAB GEOMETRY
A three-dimensional model of a solid-state electrolyte with a slab geometry containing grain and intergranular interlayer is considered. The interlayer is described by layers characterized by additional interstitial energy barriers. The interaction between mobile ions consists of long-range electrostatic repulsion and short-range attraction. The system is located between two oppositely charged electrodes. The system is simulated by the Monte Carlo method. The Coulomb energy is determined by Ewald summation with the addition "vacuum regions" on both sides of the slab. The concentration profile of mobile ions was investigated under various conditions. A comparison of the results with analytical calculations shows their quantitative discrepancy of not more than 15%.
Key words: solid-state electrolyte, grain boundary, electrode, Monte Carlo method, Ewald summation, electrical conductivity, concentration profile.
Введение. В настоящее время в электрохимических системах активно используются жидкие растворы ионных солей или полимерные ионообменные мембраны, что сопряжено с опасностью появления утечек и воспламенения. Переход к электрохимическим элементам с твердотельными электролитами может обеспечить повышение прочности, долговечности, экологичности и безопасности источников энергии, расширить диапазон рабочих температур [1-4].
Твердые электролиты являются предметом интенсивной научной деятельности ввиду широких перспектив их промышленного применения. На их основе, например, изготавливаются аккумуляторные батареи [2], топливные эле-
менты [4], суперконденсаторы [5], запоминающие устройства и т. д. Во многих случаях твердые электролиты представляют собой керамику или поликристаллиты, что требует разработки специфических методов их получения и экспериментальных исследований. Теоретические исследования таких материалов также сталкиваются с серьезными осложнениями вследствие необходимости учета дальнодействующего ку-лоновского взаимодействия, наличия межфазных границ и приэлектродных областей.
В работе [6] рассматривалась трехмерная решеточная система, в которой учитывалось только кулоновское взаимодействие. В данной статье также принимается во внимание близкодействующее ван-дер-ваальсовое притяжение.
Модель твердотельного электролита со сляб геометрией. Компьютерное моделирование реальных трехмерных кулоновских систем обычно выполняется на относительно небольшом числе частиц в пределах 100 < N < 10 000 [7, 8]. Размер системы ограничен скоростью компьютерного выполнения программы. Моделируемая система занимает весь предоставляемый ей объем. Поэтому, если в начале расчета задать координаты частиц системы в некоторой конечной области, через определенное время они способны разлететься на большие расстояния. Чтобы моделировать поведение системы при заданной плотности или давлении, необходимо поместить эти частицы в непроницаемый ящик. При этом граничные условия будут нарушать однородность системы, и для исследования макроскопических свойств потребуется значительно увеличить ее размеры, что приведет к существенному увеличению необходимых вычислительных мощностей.
Проблема поверхностных эффектов может быть преодолена путем реализации периодических граничных условий. Ящик реплицируется периодическим образом в пространстве, образуя бесконечную последовательность в трех измерениях.
В работе рассматривается трехмерная решеточная модель керамического ионного проводника размером ЬХ^ЬУ^Ь2 (моделируемая система принималось размером 30x10x10 узлов, образующих простую кубическую решетку), которая находится между двумя электродами. При этом, кроме дальнодействующего электростатического отталкивания между ионами, в модели учитывается ван-дер-ваальсовое притяжение ближайших соседей. Ионы выполняют термоактивированные прыжки в случайно выбранные вакантные узлы. Электронейтральность системы обеспечивается наличием неподвижных ионов противоположного заряда. Электроды моделируются перпендикулярными к оси х стенками, по узлам которых случайным образом разбросаны противоположно заряженные частицы со средней концентрацией ем>.
Ввиду наличия электродов возникает проблема моделирования системы со сляб геометрией, когда в двух направлениях, параллельных электродам, можно воспользоваться периодическими граничными условиями, а в третьем -система конечна. При такой геометрии для возможности использования метода Эвальда для суммирования дальнодействующих куло-новских взаимодействий с периодичностью в трехмерном пространстве необходимо добавлять «вакуумные области» с обеих сторон сляба внутри ячейки моделирования и корректировать вклады фиктивных копий системы в на-
правлении, перпендикулярном ограничивающим плоскостям сляба [9, 10]. Эффективные численные алгоритмы возможны, если коррекция производится лишь по дипольному взаимодействию, что обусловливает необходимость использования достаточно больших «вакуумных областей». Проблема заключается в корректном выборе размера «вакуумной области» и пределов суммирования в обратном пространстве. Размер «вакуумной области» при моделировании принимался равным 20 параметрам решетки.
Для воспроизводства эффекта повышенного сопротивления межзеренной границы будем моделировать ее слоем из двух плоскостей, перпендикулярных оси Х, с отличающимися межузловыми энергетическими барьерами ЦЬап-(рис. 1), где индекс I нумерует узлы решетки.
При моделировании системы случайным образом выбирается ион, находящийся в узле j. Направление его возможного перехода определяется случайным образом в один из свободных узлов ¡. Интенсивности перескоков ионов пропорциональны величине:
Wjl — ехр
(ЦЬап' - ЦЬаг, ) + J(Zj _ ^ ) + АЦСои1
квТ
(1)
где J - энергия ван-дер-ваальсового притяжения ближайших соседей; Zj, - число ближайших соседей j-го и ¡-го узла соответственно; кв - постоянная Больцмана; Т - температура.
Моделировалась одна траектория. Первые 500 Монте-Карло шагов (МКШ) отводились на установление стационарного состояния, последующее усреднение производилось по 104 МКШ.
Рис. 1. Модель энергетических барьеров вдоль оси, перпендикулярной плоскостям электродов
Разность кулоновских энергий иона до и после прыжка определяется выражением
АЦ — ттеМ _ Ц АЦСои1 - ЦСои1 ЦСои1
т^аТ:
(2)
Величину энергии кулоновских взаимодействий для системы со сляб геометрией с использованием суммирования Эвальда можно записать в виде [10]:
^Сои1 =
4Л88 0
1 N N
211
2 *=1.=1
Я.Ч,
erfc (а |г.
1 ^
4п2
пУ . а2
ехр
г .к!л
V а2 У
ео8
(к • г.)
+
+
+-
2п
V
N 2
Е у г
/=1
(3)
где у - заряд /-го иона; а = 5/Ь - параметр Эвальда; г. - радиус-вектор частицы . относительно частицы /; V = ЬХЬУЬ2 - объем системы; к = 2п(УХ / 1х , ТУ / 1у , У^ / 4 ) - вектор обратной решетки; ух, уу, у, - целые числа.
При моделировании обрезание при суммировании в прямом пространстве производилось по половине длины вдоль каждого направления, а также принималось уШаху = у'тах2 = 4, Ушах х = 12, что соответствует размеру «вакуумных областей», равному двум размерам основной ячейки моделирования в направлении оси х, и учету 334 к-векторов в обратном пространстве. При этом достигается точность около 1%, и дальнейшее увеличение области суммирования мало влияет на точность моделирования.
Результаты моделирования. С целью сопоставления результатов с аналитическими вычислениями [11] моделирование было выполнено для двух средних концентраций подвижных ионов р = 0,1 и 0,03 для системы с диэлектрической проницаемостью 8 = 41,8 и параметром решетки а = 0,4 нм при температуре Т = 1000 К.
На рис. 2 и 3 показан профиль концентрации подвижных ионов при средней концентрации р = 0,1, отсутствии/наличии притяжения ближайших соседей и отсутствии межзеренной границы.
Р
10 15 20 25 х
Рис. 2. Профиль концентрации подвижных ионов при р = 0,1, отсутствии притяжения ближайших соседей и отсутствии межзеренной границы:
1 - еМ! = 0; 2 - еМ! = 0,1; 3 - еМ! = 0,2; 4 - еМ! = 0,3
0,15
0,10
0,05
0,00
15
х
Рис. 3. Профиль концентрации подвижных ионов при р = 0,1, наличии притяжения ближайших соседей и отсутствии межзеренной границы:
1 - ^ = 0; 2 - ^ = 0,1; 3 - ^ = 0,2; 4 - ^ = 0,3
На левой положительно заряженной стенке ширина слоя пониженной концентрации растет с увеличением заряда на стенке. В то же время ионы накапливаются на противоположной стене. Слой с пониженной концентрацией появляется из-за кулоновского отталкивания от левой стенки и удовлетворения условию электронейтральности. Кроме того, учет короткодействующего притяжения между ионами приводит к более глубокому провалу концентрации вблизи правой стенки.
Такая же картина наблюдается при р = 0,03 (рис. 4 и 5). Необходимо лишь отметить увеличение ширины слоя пониженной концентрации вследствие малой объемной концентрации подвижных ионов, что требует большего объема приэлектродной области системы для обеспечения условия электронейтральности.
0 5
10 15 20 25 30
Рис. 4. Профиль концентрации подвижных ионов при р = 0,03, отсутствии притяжения ближайших соседей и отсутствии межзеренной границы:
1 - ем; = 0; 2 - ем; = 0,1; 3 - ем; = 0,2; 4 - ем; = 0,3
На рис. 6-9 показан профиль концентрации подвижных ионов при средней концентрации р = 0,1 и 0,03, отсутствии/наличии притяжения ближайших соседей и наличии межзеренной границы.
Р
5 10 15 20 25 30
х
Рис. 5. Профиль концентрации подвижных ионов при р = 0,03, наличии притяжения ближайших соседей и отсутствии межзеренной границы:
1 - ^ = 0; 2 - ^ = 0,1; 3 - ^ = 0,2; 4 - ^ = 0,3
0 5 10 15 20 25 30
х
Рис. 6. Профиль концентрации подвижных ионов при р = 0,1, отсутствии притяжения
ближайших соседей и наличии межзеренной границы:
1 - ^ = 0; 2 - ^ = 0,1; 3 - ^ = 0,2; 4 - ^ = 0,3
0,15
0,10
0,05
0,00
0 5 10 15 20 25 30
х
Рис. 7. Профиль концентрации подвижных ионов при р = 0,1, наличии притяжения
ближайших соседей и наличии межзеренной границы:
1 - ^ = 0; 2 - ^ = 0,1; 3 - ^ = 0,2; 4 - ^ = 0,3
Выбранный размер системы оказался достаточным, чтобы поведение распределения заряда вблизи межзеренной границы практически не влияло на явления вблизи стенок и наоборот.
0 5 10 15 20 25 30
х
Рис. 8. Профиль концентрации подвижных ионов при р = 0,03, отсутствии притяжения ближайших соседей и наличии межзеренной границы: 1 - ^ = 0; 2 - ^ = 0,1; 3 - ^ = 0,2; 4 - ^ = 0,3
Рис. 9. Профиль концентрации подвижных ионов при р = 0,03, наличии притяжения ближайших
соседей и наличии межзеренной границы: 1 - еМ! = 0; 2 - еМ! = 0,1; 3 - еМ! = 0,2; 4 - еМ! = 0,3
Сопоставление результатов с аналитическими вычислениями [11] показывает полное качественное совпадение результатов и количественное расхождение не более 15%, что может быть обусловлено приближениями, принятыми при аналитических вычислениях.
Заключение. Рассмотрена трехмерная модель твердотельного электролита со сляб геометрией, содержащей зерно и межзеренную прослойку, описываемую слоями, которые характеризуются дополнительными межузловыми энергетическими барьерами.
Отмечено увеличение слоя с пониженной концентрацией вблизи положительно заряженной стенки с увеличением заряда на ней. Этот эффект становится более заметным с уменьшением средней концентрации подвижных ионов.
Показано, что наличие межзеренной границы практически не влияет на распределение концентрации вблизи стенок и наоборот.
Публикация содержит результаты исследований, выполненных при грантовой поддержке научной программы Евросоюза Н0Ш20К-2020 (проект ЛМБ-734276-С0ШК) и Министерства образования Республики Беларусь.
Р
Р
Р
Р
Литература
1. Solid Oxide Fuel Cells: Materials Properties and Performance / J. Zhang [et al.] // CRC Press. 2016. 298 p.
2. High-Energy All-Solid-State Lithium Batteries with Ultralong Cycle Life / X. Yao [et al.] // Nano Lett. 2016. Vol. 16, no. 11. P. 7148-7154.
3. Ласовский Р. Н., Бокун Г. С., Вихренко В. С. Диаграммное приближение для неравновесных и неоднородных состояний решеточных систем // Труды БГТУ. Сер. VI, Физ.-мат. науки и информатика. 2010. Вып. XVII. C. 59-62.
4. Nanostructurization caused by first order phase transitions in systems with hopping dynamics / G. S. Bokun [et al.] // Solid State Ionics. 2013. Vol. 251. P. 51-54.
5. Unusual properties of a model of an intergrain boundary in solid oxide ceramic electrolytes / G. S. Bokun [et al.] // Solid State Ionics. 2017. Vol. 302. P. 25-29.
6. Ласовский Р. Н., Гапанюк Д. В., Пацаган Т. Н. Моделирование трехмерного керамического электролита с межзеренной границей // Труды БГТУ. Сер. 3, Физ.-мат. науки и информатика. 2019. № 1.C. 15-19.
7. Allen M. P., Tildesley D. J. Computer simulation of liquids. New York: Clarendon Press, 1989. 385 p.
8. Frenkel D., Smit B. Understanding Molecular Simulation. San Diego: Academic Press, 2002. 638 p.
9. Yeh C., Berkowitz M. L. Ewald summation for systems with slab geometry // J. Chem. Phys. 1999. Vol. 111. P. 3155-3162.
10. Santos A., Girotto M., Levin Y. Simulations of Coulomb systems with slab geometry using an efficient 3D Ewald summation method // J. Chem. Phys. 2016. Vol. 144. Art. 144103.
11. The effect of short-range interaction and correlations on the charge and electric field distribution in a model solid electrolyte / Т. Patsahan [et al.] // Solid State Ionics. 2019. Vol. 335. P. 156-163.
References
1. Zhang J., Hui R., Wilkinson D., Li X. Solid Oxide Fuel Cells: Materials Properties and Performance. CRC Press, 2016. 298 p.
2. Yao X., Liu D., Wang C., Long P., Peng G., Hu Y., Li H., Chen L., Xu X. High-Energy All-SolidState Lithium Batteries with Ultralong Cycle Life. Nano Lett., 2016, vol. 16, no. 11, pp. 7148-7154.
3. Lasovsky R. N., Bokun G. S., Vikhrenko V. S. Diagram approximation for nonequilibrium and in-homogeneous states of lattice systems. Trudy BGTU [Proceedings of BSTU], series VI, Physics and Mathematics. Informatics, 2010, issue XVIII, pp. 59-62 (In Russian).
4. Bokun G. S., Lasovsky R. N., Vikhrenko V. S. Nanostructurization caused by first order phase transitions in systems with hopping dynamics. Solid State Ionics, 2013, vol. 251, pp. 51-54.
5. Bokun G. S., Groda Y. G., Lasovsky R. N., Vikhrenko V. S. Unusual properties of a model of an intergrain boundary in solid oxide ceramic electrolytes. Solid State Ionics, 2017, vol. 302, pp. 25-29.
6. Lasovsky R. N., Gapanjuk D. V., Patsahan T. N. Modeling of a three-dimensional ceramic electrolyte with an intergranular border. Trudy BGTU [Proceedings of BSTU], series 3, Physics and Mathematics. Informatics, 2019, no. 1, pp. 15-19 (In Russian).
7. Allen M. P., Tildesley D. J. Computer simulation of liquids. New York, Clarendon Press, 1989. 385 p.
8. Frenkel D., Smit B. Understanding Molecular Simulation. San Diego, Academic Press, 2002. 638 p.
9. Yeh C., Berkowitz M. L. Ewald summation for systems with slab geometry. J. Chem. Phys., 1999, vol. 111, pp. 3155-3162.
10. Santos A., Girotto M., Levin Y. Simulations of Coulomb systems with slab geometry using an efficient 3D Ewald summation method. J. Chem. Phys., 2016, vol. 144, art. 144103.
11. Patsahan T., Bokun G., di Caprio D., Holovko M., Vikhrenko V. The effect of short-range interaction and correlations on the charge and electric field distribution in a model solid electrolyte. Solid State Ionics, 2019, vol. 335, pp. 156-163.
Информация об авторах
Ласовский Руслан Николаевич — кандидат физико-математических наук, доцент кафедры механики и конструирования. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: lasovsky@tut.by
Гапанюк Дмитрий Владимирович — кандидат физико-математических наук, доцент кафедры механики и конструирования. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: gapdm@mail.ru
Грода Ярослав Геннадьевич - кандидат физико-математических наук, доцент, доцент кафедры механики и конструирования. Белорусский государственный технологический университет (220006, г. Минск, ул. Свердлова, 13а, Республика Беларусь). E-mail: groda@belstu.by
Information about the authors
Lasovsky Ruslan Nikolaevich - PhD (Physics and Mathematics), Assistant Professor, the Department of Mechanics and Engineering. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: lasovsky@tut.by
Gapanjuk Dmitry Vladimirovich - PhD (Physics and Mathematics), Assistant Professor, the Department of Mechanics and Engineering. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: gapdm@mail.ru
Groda Yaroslav Gennad'yevich - PhD (Physics and Mathematics), Associate Professor, Assistant Professor, the Department of Mechanics and Engineering. Belarusian State Technological University (13a, Sverdlova str., 220006, Minsk, Republic of Belarus). E-mail: groda@belstu.by
Поступила после доработки 23.11.2019
