Содержание к диссертации
Введение
Глава 1, Обзор литературы. Описание методов компьютерного моделирования 6
I Твердые суперионные проводники: общие положения 6
II, Метод молекулярной динамики 9
III. Простые модели твердых электролитов, используемые в численном эксперименте 17
IV. Свойства и молекулярно-динамическое моделирование твердых суперионных проводников флюоритового типа: литературный обзор 21
V. Метод Монте-Карло и метод молекулярной статики 34
Глава 2. Использованные модели флюоритоподобных систем и параметры численного эксперимента 39
Глава 3 Результаты численного эксперимента: термодинамические и транспортные свойства вг^СйД^* и Bai^GdUF^ 44
Глава 4. Структурные характеристики и движение ионов фтора в системах Bai^Gd^-tt, Sri^Gd^F2+J: и SrI
Глава 5. Моделирование SrF2 и Sr^GdrF^ с применением «оболочечной модели» 59
Выводы и заключение 64
Литература 67
- Метод молекулярной динамики
- Простые модели твердых электролитов, используемые в численном эксперименте
- Свойства и молекулярно-динамическое моделирование твердых суперионных проводников флюоритового типа: литературный обзор
- Метод Монте-Карло и метод молекулярной статики
Введение к работе
Твердые электролиты (ТЭЛ), или суперионные проводники (СИП),- вещества, способные к образованию твердых фаз, где ионы одного типа находятся в разупорядоченном состоянии с высокой подвижностью и перемещаются в потенциальном поле, образованном жесткой подрешеткой ионов другого типа, -имеют важное теоретическое и прикладное значение, в том числе в квантовой электронике, прикладной оптике, сенсорной технике и т. д. В частности, в настоящее время эти материалы уже используют в малогабаритных энергоемких источниках тока, в различных датчиках составов, в солнечных элементах, таймерах и т.д.; перспективно использование твердых электролитов в устройствах прямого преобразования энергии химической связи топлива в электрическую энергию (топливных элементах).
Несмотря на интенсивный характер исследований, единой теории суперионной проводимости к настоящему времени не создано, хотя существует много различных * моделей с ограниченной областью применения. Поэтому для понимания процессов, происходящих в ТЭЛ, весьма полезны могут быть методы численного эксперимента, в частности - метод молекулярной динамики (МД), основанный (в классической версии) на численном решении ньютоновских уравнений движения для модельной системы, состоящей из сравнительно небольшого числа частиц (сотни или тысячи), взаимодействующих друг с другом согласно той или иной выбранной модели потенциала. При этом используются стандартные методы, позволяющие учесть макроскопичность реальной системы. Результаты численного эксперимента могут быть сравнены как с опытными данными, так и с результатами «строгого» теоретического исследования (например, квантовохимических расчетов). Кроме того, численный эксперимент позволяет легко и наглядно получить такие микроскопические характеристики системы (описывающие, например, пространственное распределение и движение частиц), которые в реальном эксперименте не определяются или определяются лишь опосредованно, с использованием ряда дополнительных, в той или иной мере произвольных допущений.
Целью настоящей работы является построение и изучение в численном эксперименте модели суперионной проводимости на примере фторидов со структурой типа флюорита. Методом МД моделируются стехиометрические и нестехиометрические системы флюоритовой структуры (Bab^GdJ^+j:, Sr^GdJ^-w, Sr^La^+J в широком интервале температур и составов с расчетом энергетических и транспортных характеристик моделируемой фазы и особенностей кристаллохимической структуры (характеристики пространственного расположения ионов, локальные искажения кристаллической решетки при наличии собственных и примесных дефектов). Результаты расчетов сравниваются с экспериментальными данными, что позволяет дать оценку возможностям используемых моделей.
Флюоритоподобные твердые электролиты относятся к числу СИП со сравнительно простой структурой. Поэтому модели для использования в компьютерном эксперименте для них могут быть достаточно легко разработаны и опробованы, С их помощью может быть проведено более подробное исследование процессов, происходящих в данных системах, включая те, с которыми связано практическое применение ТЭЛ. Полученные результаты, показывающие, что возможности применения простых потенциальных моделей для описания свойств исследуемых систем достаточно высоки, способствуют прогрессу исследований суперионных проводников в целом.
На защиту выносятся следующие положения.
1) При молекулярно-динамическом моделировании системы Ba^GcU^ (О < л: < 0,25) уже весьма простая модель Борна- Манера- Хаггинса позволяет получить удовлетворительное согласие между результатами расчетов и эксперимента по широкому кругу термодинамических, транспортных свойств и структурных характеристик. Подтверждено наличие суперионного перехода с «квазиплавлением» анионной подрешетки при температуре выше 1100 К. Для кристаллической фазы, существующей при более низких температурах, выявлена предпочтительная склонность ионов редкоземельной примеси и междоузельных ионов фтора к образованию тригональных дипольных комплексов по сравнению с тетрагональными. Прослежена тенденция к ассоциации примесных катионов. Воспроизведена сравнительно устойчивая простейшая кластерная структура, образующаяся в результате димеризации тригональных дипольных комплексов. В
согласии с выдвигавшимися ранее моделями, увеличение концентрации редкоземельной примеси в твердом растворе приводит к разупорядочению анионной подрешетки- локальному при малых концентрациях примеси, с тенденцией к перекрыванию разупорядоченных областей с ростом концентрации. Это перекрывание сопровождается существенным ослаблением концентрационной зависимости подвижности анионов в целом.
2) При моделировании стронцийсодержащих флюоритоподобгых систем
(Sri_xGd*F2+^ Sr^LaJ^+j) с применением аналогичной простой модели Борна-
Майера - Хаггинса также получены согласующиеся с экспериментом результаты по
предпочтительному расположению междоузельных ионов фтора, однако
вычисленные термодинамические и транспортные свойства соответствуют заметно
большей разупорядоченности анионной подрешетки в модельных системах по
сравнению с экспериментально наблюдаемой. Можно предположить, что
используемый нами модельный потенциал для взаимодействий Sr + - F" достаточно
адекватно описывает межчастичное отталкивание на малых расстояниях, но
недооценивает эффект притяжения и, как следствие, прочность кристаллической
структуры.
3) Численный эксперимент с применением «оболочечпой модели», учитывающей
поляризуемость ионов, позволяет рассчитать термодинамические и транспортные
свойства чистого фторида стронция в хорошем согласии с экспериментальными
данными. Для моделирования твердого раствора с редкоземельной примесью
(Sr^GdJ^+J требуется дальнейшее усовершенствование модельного потенциала,
поскольку имеющаяся модель для описания взаимодействий Gd+ - F^
недооценивает эффект притяжения между ионами фтора и редкоземельной примеси.
Объем и структура работы. Диссертация содержит пять глав и заключение. Глава 1 посвящена обзору литературы по экспериментальным исследованиям флюоритоподобных ТЭЛ и по применению методов численного моделирования (методы молекулярной динамики, Монте-Карло, молекулярной статики) для изучения суперионных проводников. В главе 2 описаны используемые модельные потенциалы и методики моделирования. Глава 3 посвящена изложению результатов расчетов термодинамических и транспортных свойств систем Ba^Gd^^, Sri^GdJ^+x, Sr^La^Vs с применением потенциала Борна— Майера— Хаггинса, а
Метод молекулярной динамики
В последнее время неотъемлемой частью современной науки стало численное моделирование, которое обрело широкое использование наряду с традиционными экспериментальными и теоретическими методами [16, 17]. Методы вычислительного эксперимента исходят практически из основных законов физики, что открывает дополнительные возможности в исследовании различных физических систем. Налицо также определённые аналогии при проведении физического и Затем составляется компьютерная программа, описывающая физический эксперимент и моделирующая изучаемый объект. Тестирование программы производится для проверки правильности ее работы и адекватности рассчитанных параметров с реальными физическими данными. На последнем этапе производится расчет и анализ полученных данных. Для решения подобных задач относительно твердых тел, потенциальных ТЭЛов, используется метод молекулярной динамики (МД) [18-21]. Метод молекулярной динамики (МД), наряду с методом Монте-Карло (МК), принадлежит к методам молекулярного моделирования, проводимого с использованием современных компьютеров, которые с успехом применяются при исследовании конденсированного состояния вещества. Преимуществами этих методов является возможность одновременного получения как детальнейшей информации о структуре изучаемых систем, так и их термодинамических характеристик. Часто методы молекулярного моделирования называют численным или компьютерным экспериментом в статистической и молекулярной физике. Действительно, результатом применения данных методов всегда являются конкретные численные значения исследуемого свойства системы при заданных значениях внешних параметров-
Кроме того, с точки зрения чисто модельных систем, таких, как твердые сферы, решеточные модели и др_, рассматриваемых в статистической физике, численный эксперимент является единственным источником достоверных данных для проверки всевозможных теоретических приближений. Особенно актуальным применение методов компьютерного моделирования становится при изучении объектов или состояний, малодоступных реальному физико-химическому эксперименту. Поэтому эти методы занимают промежуточное положение между реальным экспериментом и теорией. Впервые этот подход был применен Олдером к системе твердых дисков [22]. А в работах Рамана [23] и Верле [24] уже рассматривались более реалистические системы с потенциалом Леннард-Джонса. Подробное изложение методов молекулярной динамики с обзором основных достижений в исследовании свойств конденсированных фаз можно найти в литературе [16, 25-27]. Поскольку в настоящее время имеется достаточно большое число различных модификаций метода молекулярной динамики, приспособленных для изучения разнообразных систем в разных условиях, то говорят уже не об одном методе, а о целой группе методов молекулярной динамики, объединенных основной идеей описания системы на основе численного решения уравнений движения молекул или атомов, составляющих данную систему, при заданных начальных и граничных условиях,
В результате численного интегрирования системы дифференциальных уравнений получают траекторию системы в фазовом пространстве, т.е. совокупность координат гі, Г2, Гз, ... Гдг и скоростей Vi, v2, V3J .-. vw всех N молекул в зависимости от времени f. Это позволяет в рамках метода МД описывать не только равновесные свойства, но также характеризовать релаксационные процессы и получать информацию о кинетических коэффициентах (например, диффузии). Как известно из статистической физики, полный набор координат и скоростей молекул системы в данный момент времени характеризует микроскопическое состояние системы. Средние значения термодинамических свойств системы получают путем усреднения но времени, или по фазовой траектории, соответствующей величины, зависящей от микроскопического состояния. Например, для некоей величины A (r f), г2(Д ґз№ -- Ы Х vi( ) v2(0 v3(0, VAO) і) среднее значение А, измеряемое в эксперименте, равно Іде для краткости записи введено следующее обозначение для совокупности векторов, характеризующих микросостояние системы молекул:
В методе МД фазовую траекторию и, следовательно, любую динамическую величину A(TX(t)t г2(0, г3(/) ...rN(0, vi(f), v2(r), v3(0t - vN(0 0 рассчитывают, численно интегрируя юіассические уравнения движения молекул. Естественно, что при таком подходе можно иметь дело только с конечным отрезком фазовой траектории, точнее, с ее приближенным отображением. Поэтому в методе МД необходимо, с одной стороны, рассчитывать отрезок фазовой траектории достаточной протяженности, чтобы получить представительные результаты измерения. С другой стороны, этот расчет должен быть по возможности более точным, чтобы утверждать, что усреднение по времени соответствует усреднению именно по фазовой траектории. Очевидно, что чем точнее требуемый расчет, тем больше времени будет затрачено на получение отрезка траектории необходимой длины. Удовлетворение этим условиям требует некого компромиссного решения, что нужно иметь в виду при выборе конкретного алгоритма для численного интегрирования уравнений движения. При численном решении уравнений движения мы будем иметь набор точек на данном отрезке фазовой траектории. Если точки разделены равным промежутком времени 5?, то среднее значение величины A({r(f) v( )} 0 по отрезку траектории с п точками равно
Простые модели твердых электролитов, используемые в численном эксперименте
Многие ионные и суперионные проводники могут рассматриваться как кристаллы с чисто ионным характером связи, что позволяет использовать простую модель «жестких ионов» [31-33]. При этом потенциальная энергия решетки аппроксимируется суммой парных взаимодействий между ионами. Потенциалы ионного взаимодействия определяются различными методами. Межатомный потенциал может рассчитываться из первых принципов на основе квантовомеханических представлений в рамках модели электронного газа [34], или с помощью полуэмпирических квантово-механических уравнений для модели Малликена [35]. Получаемые результаты либо используются в табличном виде, либо с их помощью находятся параметры модельных уравнений, описывающих межчастичные потенциалы. В рамках модели жестких ионов (rigid-ion) потенциал Борна -Майера -Хаггинса в общем виде записывается следующим образом [36]: эффективные заряды взаимодействующих ионов І nj, находящихся на расстоянии г, е - заряд электрона, гц9 nj -числа электронов на внешних оболочках, ОЇ, GJ - эффективные ионные радиусы, pi, pj -электронные поляризуемости ионов, Ьф pij и Су -постоянные.
Первое слагаемое в (1) описывает кулоновское взаимодействие ионов, второе - энергию отталкивания за счет перекрывания электронных оболочек, третье - энергию поляризации, четвертое - энергию притяжения ван-дер-ваальсовского взаимодействия. В рамках модели жестких сфер используется еще один подход, в котором рассматриваются «сферы с мягкой сердцевиной» (ionic soft-core system) [14,36]. При этом для описания межчастичного взаимодействия используется потенциал Митры: В рассматриваемом подходе межчастичный потенциал состоит из двух составляющих, характеризующих электростатическое (кулоновское) взаимодействие и «отталкивание» сердцевин. Потенциал Борна - Майера - Хаггинса часто вводится в упрощенном виде: где Ац - постоянные. Второе слагаемое в (24) описывает энергию отталкивания за счет перекрывания электронных оболочек. Учет поляризации в компьютерном моделировании требует увеличения времени расчетов, однако он необходим при вычислении энергий дефектов в полярных кристаллах. В ранних работах по этой теме рассматривались модели точечных диполей, в которых динольный момент pi для точечного диполя считался пропорциональным величине локального напряжения электростатического поля: /а=аЕ; такая модель приемлема для изолированных или небольших молекул, но неприменима при расчетах кривых дисперсии фононов для ионных кристаллов и (статических) диэлектрических свойств [37,38], Такие погрешности возникают в результате пренебрежения тем фактом, что энергия короткодействующего отталкивания связана с поляризацией, jha связь становится очевидной, если рассматривать физический источник поляризации атома как искажение его элекгронного облака, т.е. использовать оболочечную модель (Дика- Оверхаузера [39]), схематически показанную на рис. 3. Согласно модели поляризуемые электроны с валентных орбиталей представляются невесомой или обладающей малой массой оболочкой, которая связана «пружиной» с ядром; в ядре сконцентрирована вся масса атома или ее большая часть. Поляризация описывается как смещение оболочки относительно ядра, сопровождающееся растяжением «пружины» (энергия растяжения рассчитывается по закону Гука); так как короткодействующее отталкивание действует между оболочками, то модель включает в себя требуемую связь между отталкиванием и поляризацией. Оболочечная модель позволила успешно вычислять механические (пластические), диэлектрические и динамические свойства решетки твердых тел. Общий заряд катиона От + щ Общий заряд аниона Q2+ ( Заряд обояочш ( Заряд оболочки О У ие. J. Схематичное представление о&олочечной модели для катионов и анионов в ионном кристалле. qJf q2 - заряды ядер; Qj, Q2 -заряды оболочек; R - расстояние между ионами; W]t W2 -величины поляризации. Все модели с использованием вышеописанных потенциалов включают в себя множество подгоночных величин, таких как параметры короткодействующего отталкивания, поляризации и пр. Они могут быть определены либо из известных физических свойств (эмпирическая процедура), либо с помощью теоретических моделей. Оба подхода широко обсуждаются в литературе [40-42], Принципиальным ограничением метода эмпирической подгонки параметров является то, что, например, для многих «сложных» материалов точно определена лишь кристаллическая структура.
Другие физические характеристики неизвестны, поскольку для их измерения требуется наличие монокристаллов высокого качества. В некоторых случаях трудно использовать даже структурные данные, в частности, для нестехиометрических или примесных материалов [43]. Кроме того, в экспериментальном определении значения энергий решетки большой вклад вносят взаимодействия, не учитываемые моделью, поэтому такие данные также должны быть применены с осторожностью. Подобное обсуждение проблем эмпирического подбора потенциалов взаимодействия можно найти в [44]. В ряде работ [45-47] была показана возможность вычисления потенциалов короткодействующего отталкивания в рамках модели электронного газа. Вероятно, большая точность может быть достигнута путем использования метода Хартри Тшпия; жсперижішльїшс данино для ютгорьж [54 57] указывают т то. что при t&mMpm урс на 200 фазовый переход, соп хмжтрешровпдиосш кристалла 1\л& указанных веществ в твердом coca исследования транспортных и структурных характеристик, включая прямые измерения электропроводности [56,58-61]. В [61] специально изучалась электропроводность при высоких давлениях, в частности, для /?-PbF2 - кубической флюоритоподобной модификации фторида свинца, нахождение коэффициента диффузии анионов с использованием методов ЯМР [60,62], анализ транспортных и других свойств с применением метода термически стимулируемой деполяризации (для чистых флюоритоподобных фторидов и твердых растворов, содержащих также трехвалентные катионы) [63,64], нсйтронографические исследования [65,66]; в [66] данные по дифракции нейтронов применялись для изучения температурной зависимости параметра ячейки кристалла и теоретического расчета величины теплоемкости /?-PbF2, которая затем сравнивалась с данными калориметрических измерений Согласно данным, приведенным в [66], теплоемкость при постоянном давлении Ср для /?-PbF2 составляет при температурах 500-550 К (соответствующих кристаллической фазе без анионного разупорядочения) примерно 10R (R-универсальная газовая постоянная), при 800 К (когда стабильна фаза с анионным разупорядочением)- около ПК, а в области суперионного перехода (711 К) достигает 25-ЗОД (теплоемкость при постоянном объеме Су может быть получена путем вычитания из Ср составляющей, связанной с тепловым расширением кристалла» которая по расчетам авторов [66] равна примерно 0,8Д для твердой фазы без разупорядочения анионной подрешетки, около R- для суперионной фазы, а в области суперионного перехода доходит до 2R). Важной особенностью нестехиометрических фторидных материалов служит неизменность их кристаллической струиуры при варьировании состава. Конечно, при этом образуются структурные дефекты, которые могут приводить даже к возникновению ближнего порядка в некоторых структурных единицах (захватывая область до нескольких постоянных решетки), однако новой, полностью упорядоченной фазы не возникает. Следует отметить, что при большой концентрации дефектов они часто образуют агрегаты (из точечных дефектов), называемые кластерами. На основе нейтроио- и рентгенографических исследований нестехиометрических флюоритовых фаз M RJ -kt (М=Ва, Si\ Са, Cd, R= La-Lu, Y) были предложены две модели, описывающие дефектную структуру кристаллов - двухпозиционная [68,69]
Свойства и молекулярно-динамическое моделирование твердых суперионных проводников флюоритового типа: литературный обзор
В ряде работ [45-47] была показана возможность вычисления потенциалов короткодействующего отталкивания в рамках модели электронного газа. Вероятно, большая точность может быть достигнута путем использования метода Хартри Тшпия; жсперижішльїшс данино для ютгорьж [54 57] указывают т то. что при t&mMpm урс на 200 фазовый переход, соп хмжтрешровпдиосш кристалла 1\л& указанных веществ в твердом coca исследования транспортных и структурных характеристик, включая прямые измерения электропроводности [56,58-61]. В [61] специально изучалась электропроводность при высоких давлениях, в частности, для /?-PbF2 - кубической флюоритоподобной модификации фторида свинца, нахождение коэффициента диффузии анионов с использованием методов ЯМР [60,62], анализ транспортных и других свойств с применением метода термически стимулируемой деполяризации (для чистых флюоритоподобных фторидов и твердых растворов, содержащих также трехвалентные катионы) [63,64], нсйтронографические исследования [65,66]; в [66] данные по дифракции нейтронов применялись для изучения температурной зависимости параметра ячейки кристалла и теоретического расчета величины теплоемкости /?-PbF2, которая затем сравнивалась с данными калориметрических измерений Согласно данным, приведенным в [66], теплоемкость при постоянном давлении Ср для /?-PbF2 составляет при температурах 500-550 К (соответствующих кристаллической фазе без анионного разупорядочения) примерно 10R (R-универсальная газовая постоянная), при 800 К (когда стабильна фаза с анионным разупорядочением)- около ПК, а в области суперионного перехода (711 К) достигает 25-ЗОД (теплоемкость при постоянном объеме Су может быть получена путем вычитания из Ср составляющей, связанной с тепловым расширением кристалла» которая по расчетам авторов [66] равна примерно 0,8Д для твердой фазы без разупорядочения анионной подрешетки, около R- для суперионной фазы, а в области суперионного перехода доходит до 2R).
Важной особенностью нестехиометрических фторидных материалов служит неизменность их кристаллической струиуры при варьировании состава. Конечно, при этом образуются структурные дефекты, которые могут приводить даже к возникновению ближнего порядка в некоторых структурных единицах (захватывая область до нескольких постоянных решетки), однако новой, полностью упорядоченной фазы не возникает. Следует отметить, что при большой концентрации дефектов они часто образуют агрегаты (из точечных дефектов), называемые кластерами. На основе нейтроио- и рентгенографических исследований нестехиометрических флюоритовых фаз M RJ -kt (М=Ва, Si\ Са, Cd, R= La-Lu, Y) были предложены две модели, описывающие дефектную структуру кристаллов - двухпозиционная [68,69] и однопозиционная [67,70-72]. Согласно первой модели дополнительные ионы фтора располагаются одновременно на осях симметрии второго и третьего порядка, а во второй модели - на осях только второго или только третьего порядка. Следует отметить, что данные разных авторов не всегда согласуются между собой (см., в частности, краткий обзор в [71]). Тип и природа кластеров, минимальная концентрация примеси для начала кластерообразования до настоящею времени также трактуются неоднозначно. Так, например, согласно Хартогу [73,74], кластеры образуются только в твердых растворах Sr R -w, гДе R-Er-Lu,Y, что несколько расходится с практикой. Если межузельные ионы фтора располагаются на тройных осях симметрии (позиции 32/), то в твердых растворах образуются кластеры. Процессы образования различных кластеров во флюоритовых матрицах являлись предметом детальных расчетов, проведенных Кэтлоу с соавторами [75-79]. Для объяснения особенностей поведения концентрационных зависимостей характеристических параметров твердых растворов MI RJ -K Ивановым-Шиц с соавторами была предложена модель дефектных областей [80], развивающая представления перколяционной модели Хартога и модели «ускоренного ионного переноса», предложенной Скунманом [81]. При изоморфном введении добавки RF3 в кристаллическую матрицу MF2 образуются «дефектные области» (домены), содержащие структурные дефекты кристаллической решетки флюорита: редкоземельные ионы и компенсирующие их заряд ионы фтора. Дефектная область имеет состав и локальную кристаллическую структуру, отличающиеся от исходной матрицы. За пределами «дефектной области» как структурные искажения флюоритовой матрицы, так и энергетические возмущения, вызванные ядром, сведены к минимуму. Иными словами, можно считать, что нестехиометрические флюоритоподобные твердые растворы состоят как бы из двух «фаз» - матрицы MF2 и «дефектной фазы». Возможность существования двух «фаз» в концентрированных твердых растворах подтверждается данными ЯМР [82-84] и спектроскопических оптических исследований [85,86], которые указывают на образование в кристаллической решетке как областей с повышенным содержанием примесных катионов, так и областей, по структуре близких к чистой флюоритовой матрице.
Предметом пристального исследования стали электролитические свойства твердых растворов типа Ba R +x, для которых характерно образование кубооктаэдрических кластеров АбР3б-з8 (где А- Ва или R). В рамках модели дефектных областей будем считать, что ядром дефектной области служит кубооктаэдрический кластер, а дополнительные ионы фтора, расположенные вне кластера (Fcxt), находятся в периферии дефектной области и занимают позиции внутри пустых анионных кубов флюоритового типа твердого раствора. В зависимости от концентрации Fejt будет меняться проводимость дефектной области, а, следовательно, и всего кристалла в целом. При интерпретации экспериментальных данных для реальных физико-химических систем (в том числе и для твердых электролитов) весьма интересен вопрос о степени применимости относительно простых «классических» моделей межчастичного взаимодействия для описания их поведения, В решении этого вопроса большую роль могут сыграть методы численного эксперимента- метод Монте-Карло (МК) и метод молекулярной динамики (МД), которые позволяют, исходя из достаточно простых модельных потенциалов, рассчитать широкий набор термодинамических и структурных свойств моделируемых систем. Метод МД также дает возможность получить транспортные характеристики и воспроизвести эволюцию системы во времени. В численном эксперименте могут быть рассчитаны такие характеристики поведения системы на уровне атомов, молекул или ионов, прямое экспериментальное определение которых либо затруднено, либо даже невозможно. Результаты «численного эксперимента» могут быть сравнены как с реальным экспериментом, так и с данными сравнительно строгих по своим исходным посылкам теоретических подходов (использующих, например, квантовые расчеты ab initio). Краткое описание метода МД было дано в II. Он использовался при изучении широкого круга СИП :AgCl, Ag3SI, L13N, CaF2 ,ZrC 2 и ряда других [35,87,88]. Необходимо специально отмстить работы [21,89,90], посвященные кристаллам с флюоритовой структурой. Так, в [89] в молекулярно-динамическом численном эксперименте весьма подробно изучен характер анионного транспорта в кристаллах CaF2, SrF2, PbF2 при температурах как ниже, так и выше точки суперионного перехода; был сделан вывод о том, что в разупорядоченной анионной подрешетке
Метод Монте-Карло и метод молекулярной статики
Параллельно с методом МД в расчетах различных состояний ТЭЛ используется также метод Монте-Карло (МК), суть которого заключается в следующем. Состояние любой системы характеризуется значениями числа частиц N входящих в систему, объемом V и полной энергией Е этой системы. На микроскопическом уровне существует множество различных способов (или конфигураций), в которых может реализовываться данное микросостояние (JV,V,E). Если рассматриваемая система находится в контакте с тепловым резервуаром с заданной температурой Г, то из курса статистической физики известно, что вероятность Ps того, что система находится в микросостоянии s с энергией Es описывается выражением где суммирование происходит по всем Ммикросостояниям системы. Среднее (по ансамблю) значение физической величины А дается соотношением где т - ограниченное число из полного числа М конфигураций системы. Для того чтобы смоделировать бесконечно большую физическую систему, будем рассматривать бокс с периодическими граничными условиями, как это сделано в методе МД. Вычислительная процедура метода Монте-Карло (МК) состоит из нескольких этапов [108,109]. Произвольным образом с помощью генератора случайных чисел выбирается частица, и затем ее координаты г изменяется в соответствии с уравнением где R- тройка случайных чисел из интервала (0,1) и S- максимально допустимое смещение частицы (порядка 0,01 нм).
После этого рассчитывается конфигурационная (потенциальная энергия) и полученное значение Е сравнивается с величиной энергии предыдущей (стартовой) конфигурации Е. Если АЕ=Е -Е 0, то определяются необходимые физические параметры и осуществляется переход к новой конфигурации. Если Д 0 , то вычисляется вероятность перехода =схр(-Е/кТ). Величина сравнивается со случайным числом R1 (из интервала (0,1)) и в случае Rf новая конфигурация принимается. Если Я\ то система возвращается к предыдущему состоянию. Такая процедура, называемая методом выборки Метрополией повторяется достаточное число раз, и затем вычисляются средние по конфигурациям, которые статистически независимы друг от друга. В силу эргодической гипотезы следует, что результаты усреднения метода МК и усреднения по времени метода МД одинаковы. Методом МК были изучены многие ТЭЛ. Ниже рассмотрим результаты расчетов для флюорита CaF2 [ПО], т.е. соединения, изоструктурного и сходного по своим свойствам с BaF2. На рис, 8 показана гистограмма реального распределения ионов около центра гранецентрированной кубической решетки при двух температурах. В случае идеальной флюоритовой структуры все анионы находятся на расстоянии 2,55 А от центра куба. При 1773 К анионы симметрично распределены в окрестностях узлов, но при 1553 К они в большей степени смещены к междоузельным позициям. Рис. 8. Функция распределения анионов фтора относительно центра куба катионной подрешетки при Т=1773 ЦІ), Т=1553 К (2), Теплоемкость кристалла была вычислена в соответствии с соотношением Как видно из рис. 9, в CaF2 существует фазовый переход при 1470 К, что качественно согласуется с реальным физическим экспериментом. Для определения коэффициента диффузии необходимо знать временные зависимости среднеквадратичных смещений, однако в МК-вычислениях можно определить «эффективный коэффициент диффузии». На рис. 10 показаны временные зависимости среднеквадратичных смещений для катионов и анионов. Из него видно, что анионы участвуют в диффузионном движении, в то время как катионы колеблются около узлов решетки.
Таким образом, основные характеристики CaF , как и следовало ожидать, аналогичны характеристикам BaF2. Описанные выше методы МД и МК требуют больших затрат машинного времени. Более простыми и в общем случае более экономичными по отношению к компьютерным ресурсам являются методы статического моделирования, применимые к широкому спектру различных систем [40]. В результате статического моделирования может быть получена величина энергии системы в любой точке фазового пространства системы, а это позволяет находить такие важные характеристики системы, как пространственную конфигурацию стабильных и метастабильных состояний систем, разность энергий системы в этом состоянии. Иными словами, удается рассчитать как энергию образования идеальной кристаллической решетки, так и энергию кристалла с дефектами. Вычисление энергии системы производится при нуле градусов Кельвина, тем не менее, как было показано Гилланом [lllj, такне расчеты являются хорошим приближением энтальпий образования кристаллов из свободных атомов и при более высоких температурах. При возникновении в твердом теле дефекта, вокруг него происходит релаксация решеточных ионов, и в процессе вычислений ионы занимают позиции, отвечающие минимуму потенциальной энергии. Энергии релаксации достаточно велики в ионных кристаллах ввиду сильного взаимодействия заряженных дефектов (или мигрирующих ионов) с решеткой. Методы расчета характеристик твердых тел были предложены Моттом и Литглтоном [112] и развиты Лидьярдом и Норгетом [113].
Согласно этому подходу считается, что кристалл можно представить в виде двух областей; внутренней (1) и внешней (2) . Внутренняя область содержит дефект и ближайшие к нему атомы (100-500 частиц), ее структура и энергия рассчитываются точно, используя потенциалы взаимодействия ионов. Общая энергия может быть представлена [114] в виде суммы трех слагаемых ь Е2, Е связанных с энергиями области 1, взаимодействия области 1 и 2, и области 2: Здесь х и а - координаты атомов во внутренней и внешней областях соответственно. Энергии Е2 и Е$ могут быть вычислены с использованием приближенных методов. Используются, в частности, методы градиентного спуска; из них весьма эффективен мегод сопряженных градиентов [115]. Есть и другие процедуры минимизации [ 116,117] Работы по моделированию структур различных минералов, осуществленные Паркером и сотр. [118], Сандерсом и сотр. [119], представляют собой хорошую иллюстрацию применения метода статического моделирования к системам с большим количеством переменных. Наиболее часто употребляется потенциал центрального двухцентрового взаимодействия, хотя ограничения такой модели становятся все более очевидными [ 120-122]: введение трех- или более высокоуровневых многоцентровых взаимодействий в коды компьютерного моделирования вызывает серьезные трудности, В отличие от методов МД и МК, в методах молекулярной статики при рассмотрении межатомных потенциалов, кроме электростатической [123-127] и короткодействующей [128-130,37] составляющих, часто учитывается составляющая, отвечающая за атомную поляризацию (см. гл. 2).
