Спросить
Войти
Категория: Математика

ПРИМЕНЕНИЕ СЕТОЧНЫХ МЕТОДОВ РЕШЕНИЯ ЭВОЛЮЦИОННЫХ ЗАДАЧ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ИНДУЦИРОВАННОЙ ЭЛЕКТРОННЫМ ЗОНДОМ ЗАРЯДКИ ДИЭЛЕКТРИКОВ

Автор: Масловская Анна Геннадьевна

УДК 519.6: 537.9

А.Г. Масловская, А.А. Красновид, А.В. Сивунов

ПРИМЕНЕНИЕ СЕТОЧНЫХ МЕТОДОВ РЕШЕНИЯ ЭВОЛЮЦИОННЫХ ЗАДАЧ ДЛЯ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ИНДУЦИРОВАННОЙ ЭЛЕКТРОННЫМ ЗОНДОМ ЗАРЯДКИ ДИЭЛЕКТРИКОВ

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

The paper presents the results of dynamic process simulation caused by charging-up effect in polar dielectrics under the investigation with the scanning electron microscope. The calculation is based on the joint solving of the continuity equation and Poisson equation taking into account the intrinsic radiation-induced diffusion process in irradiated specimen. The initial electron distribution is determined by Monte-Carlo method. The computational is performed with programming in Matlab application package using finite deference and finite element methods. The mathematical model of the nonstationary charging process in polar dielectrics allows the charge density and potential distribution to be estimated at the given experimental parameters.

Введение

В настоящее время методики растровой электронной микроскопии (РЭМ) активно применяются в качестве аналитических методов исследования и технологической обработки различных полярных материалов. Это обусловлено широким спектром эффектов, возникающих при воздействии электронных пучков на полярные диэлектрики, а также высокой чувствительностью таких материалов к электрическим и тепловым воздействиям электронного зонда, что позволяет получать отклик и создавать способы формирования изображения и исследования электрических свойств образцов [1].

При электронной бомбардировке диэлектрического образца на его поверхности и в объеме поглощенные электроны могут образовывать локальные заряженные области, способные нерегулярным образом отклонять первичный пучок. Эффект зарядки возникает при любых увеличениях РЭМ и любых реальных токах зонда. Исследованию вопроса зарядки диэлектрических материалов при изучении методами РЭМ посвящен ряд работ как теоретического, так и экспериментального характера [2-7]. При этом все чаще в практике подобных исследований применяют средства и методы математического моделирования. Так, в [45] предложено использовать для определения объемного распределения зарядов в образце метод Монте-Карло. Активно применяются и детерминированные модели, описывающие явление зарядки диэлектриков [5-7]. Однако многие аспекты изучения процессов накопления и релаксации зарядов в диэлектрических материалах остаются неизученными. Цель данной работы - разработка, построение схемы численного решения и программная реализация математической модели, описывающей эффект зарядки, возникающий в процессе облучения полярных диэлектриков электронными пучками средних энергий.

Физико-математическая модель

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

где р - объемная плотность распределения зарядов, Кл/м3; ] - плотность тока, А/м3; а - удельная электрическая проводимость, 1/(Омм); Е - поле, В/м; е - удельная диэлектрическая проницаемость среды, 1; е0 - диэлектрическая постоянная, Кл/(В-м); в - генерационное слагаемое.

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

= -^гайр,

где П = /п— - коэффициент диффузии электронов, м2/с; /п = VГ / Е = е • I /(т • V) - дрейфовая

подвижность электронов, м2/(Вс); Е - напряженность электрического поля, В/м; Vаг - дрейфовая

скорость электрона, м/с; I - средняя длина свободного пробега, м; V - средняя тепловая скорость, м/с; т - масса электрона, кг; к - константа Больцмана, Дж/К; Т - температура, К; е - заряд электрона, Кл.

Если облучение поддерживается достаточно длительное время, то электроны, участвующие в диффузионном токе, создают объемные заряды, поле которых противодействует диффузии. Дрейфовую плотность тока электронов можно выразить через собственную радиационно-стимулированную проводимость а = /& Р и электрическое поле Е,

создаваемое зарядом: \\(г =аЕ. Таким образом, плотность тока проводимости представляется суммой диффузионного и дрейфового токов \\ = \\+ \\ (г .

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

Поскольку уравнение Пуассона носит локально-мгновенный характер, с его помощью можно определить распределение потенциала ф в зависимости от мгновенных положений зарядов. Принимая во внимание цилиндрическую симметрию задачи и учитывая, что (щ = айпЕ + Egradа - ПАр,

Рис. 1. Геометрия модельного образца.

итоговая система, формализующая математическую постановку задачи, описывается совокупностью соотношений

р ^р2-т• íЕдр+е2р+в

2

д 2р 1 др д 2р

2 +--+ 2

дг г дг дг

д 2 ф 1 дф д 2 ф р

-ч 2 + = - , Е = -§гаФ

дг г дг дг ££0

где 0 £ г £ Я,0 £ г £ г,0 £ г £ Т .

Для замыкания математической формулировки диффузионное уравнение и уравнение Пуассона, выраженные системой (2), необходимо дополнить начальным

р(г, г, го )=р0 (г, г) (3)

и, соответственно, граничными условиями:

^ = о, р = 0, ф

дг _=п & ^=г,г =Я & дг

= о, ф

2=х ,г = Я

Схема численного решения задачи

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

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

= {г = гЛ, 1 = 0, п; г- = ]Н2, ] = 0, т; гк = кт,т = 0,1,2,...},

на этой сетке будем аппроксимировать дифференциальную задачу, причем на шаг по времени накладываем ограничение, определяемое условием Куранта [9]: 1

2(^12+%) ■

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

т£12 +г

-т р)2-т Е),

р+1,] р-1,]

2И1

+ в

р,+1,, - р, +р,-1.

кк L +1 р+1, - -р&~

2к1

р(г,, г], г0 ) = рр, ] , р] = 0, - Зр] + 4р^12

к+12 р2, ]

2к1

для второго дробного шага имеем

1 =рк+12 +г
1/2 -т Е),,

рк+1/2 - к+п пк+\\/2 - 2 Пк+У2 + рк+^2 Л,,+1 Л,у-1 , г. А7/,у+1 2г /,, ^ Н/,у-1

Р(г,, г,, (у 2 )=Р1

+ О

гт = о,

ЗрРО1 + 4р£1 -Ра

0.

Но чтобы рассчитать распределение плотности заряда, необходимо найти значение напряженности в узлах расчетной области по формуле Е = - grad ф.

Решение уравнения Пуассона для интересующего момента позволяет рассчитать поверхность распределения потенциала по объему исследуемого образца:

22

Э2ф + 1 Эф + Э ф = р

Эг2 г Эг Эг2 ££с

г=0, г=0

ф1е "1 г1 г1 а1

ф2 > = 1 г2 г2 а2

фзе 1 г3 гз _ аз

Ф = о.

Мг=Я, г=2

Реализация последней задачи проведена на основе метода конечных элементов.

Разобьем расчетную область ^ на подобласти в виде треугольников, как показано на рис. 2. Пронумеровав все узлы, получим разбиение всей области на Е = 2 х( N - 1)х(М -1) треугольных конечных элементов с общим числом нумерованных узлов Р = N х М, где N - количество узлов вдоль оси г, а М -количество узлов вдоль оси г. Обозначим границу, на которой заданы условия Дирихле, Гф = Г1 + Г4, а на границе с условиями

Неймана - Гд = Г2 + Г3.

Согласно [10] интерполяционная формула для некоторой рис 2. Разбиение области функции ф на треугольном линейном элементе е, используя на треугольные элементы. значения переменных г и г в узлах элемента, приобретает вид

"1 г 1 г 1г

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

показано на рис. 3. Индексы / = 1,3 указывают на Рис. 3. Линейный треугольный элемент. порядок обхода узлов в элементе.

Вычисляя обратную матрицу, получим так называемую функцию элемента:

фе = не (г, г ф + не (г, г )ф2е + не (г, г )ф

2
2
7.
1
2
3
3
1 л

а = —

2
1 Г
1 Г ¿2 1 г3 г3

площадь треугольного элемента;

Н1 = 2Ае [(г2 ¿3 - Г3¿2 )+ (г2 - ¿3 )г + (г3 - Г2 У];

Н = 2Ае [(г3г1- г23)+(г3- г)г +(г - г]; Н = 2А ^- гг)+^- ¿2)г+(г- г] •

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

нумерованных узлах, т. е. в форме

9(Г, ^ )»j(r, 7 ) = ^фрНР (r, * ), (6)

где Р - число нумерованных узлов.

Система базисных функций Нр (г, z), р = 1, Р обладает свойством полноты, поскольку при

Р—ют решение может сколь угодно точно аппроксимировать искомую функцию. Далее, подставив приближенное решение (6) в дифференциальную задачу (5), получим не тождественный нуль, а некоторую функциональную невязку Яп(г,z) по расчетной области О и невязку Лг (г,z) - на

границах Г2 и Г3.

п / ч Э2ф 1 Эф Э2 ф р(г,z) п / ч Эф Яп(г,2) = —^ + —+ + , (г,^) = ^т . Эг г Эг Эz ее0 " Эп

В соответствии с методами взвешенных невязок требуем ортогональности этих невязок и специальным образом подобранных весовых функций Щ.(г,z), 5 = 1,Р для невязки ЯО(г,z) и Щ,(г,z), 5 = 1,Р - для невязки Лг (г,z). Для непрерывных функций ЯО(г,z) и Лг (г,z) это означает равенство нулю скалярных произведений или Яп(г,z)+Лг (г,z) = 0, что приводит к

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

Щ )+(«г, Щ )={( ф + Г ЭФ + ф + )-Щ. (г •& + Д Эф ) Щ (г,, >Г = 0,, = .

где Щ, (г, z) - весовые функции для внутренних узлов расчетной области; Щ, (г, z) - весовые

функции для граничных узлов расчетной области; 5 = 1, Р . Весовые функции выбираются таким образом, чтобы они были ортогональны невязкам по области по границе. В соответствии с методом Галеркина взвешенных невязок весовые функции равны базисным: Щ = Н. ,Щ = Н..

Ослабим формулировку задачи, используя первую формулу Грина. Весовые функции Н. на границе Гч можно выбрать таким образом, чтобы интеграл по границе Гч сократился с частью интеграла по границе Г, то есть принять Н. = -Н. . Тогда получим выражение:

I ,1 - + р=1 р I О V Эг Эг Эz Эz

Уфр |-|

=-,Г р(г, z)

Н Н р Н Н р )

¿О + I Н •/

1 ЭНр)

О V ее0

¿О- 1(Н. -^¿г,. = 1, Р.

у гЧ . Эп)

Неоднородная система линейных алгебраических уравнений порядка Р (7) в векторно-матричной форме имеет вид

Кф = /, (8)

где элементы к.р матрицы К и правых частей / образованы суммированием вкладов отдельных конечных элементов:

к = у ке А = у Г

е=1 е=1

СЛАУ (7) или (8) называют глобальной СЛАУ, матрицу К и вектор правых частей /соответственно глобальной матрицей (или матрицей жесткости) и глобальным вектором правых частей.

Для отдельного конечного элемента Ое +ге с локальными номерами узлов 1, 2, 3 на основе глобальной СЛАУ, элементы матрицы и вектора правых частей которой определяются выражением (8), можно построить локальную СЛАУ Кефе = /е. В матричном виде получим

" н )

Эг ЭН3е

ЭН7 ЭН2е ЭН3е)

Эг Эг

( эн7 )

Эz ЭН2

дне эн 2 эн

Эz Эz Эz

1 г

( н7) не

Ч Н3 у

( Н ЭН2 дН3)

Эг Эг Эг

Ре (Г z)

Ое ££0

Эф Эп

¿г, е = 1, Е.

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

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

,е _ | 1 ^2 - ¿3 ) ( ((г2¿3 - Г3¿2 )+(z2 - ¿3 )г + (г3 - Г2 )z)"

1 = I

г 2 Ае

2 Ае

(z2 - ¿3 )2 +(г3 - Г2 )2 )

(2 А7 )2

С учетом того, что |¿О = Ае и в соответствии с теоремой о среднем

= 1 (г1 + Г2 + Г3 ) , ^ = 1 (zl + ¿2 + ¿3 ) ^

получаем кЦ - —I —-3 I—--2-3 ^ e 3—-— . По аналогии находим остальные коэффициенты

1Г - ) 1 _ (г2 _ г3 )2 +(гз - г2 )2 3Г 2 ) г* 4

матрицы жесткости.

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

Ре (r, z ) ( ((r2 Z3 - r3 z2 )+(z2 - Z3 У + (r3 - r2 )z) \\jW - 1 Р

- Г^-^ I w^3 -3--/ ■ v-- -ъг ■ v-з -ir, ^dW---^-■ Ae J ee0 ^ 2Ae j 3 ee0

где pr =pe (i, z) .

Тогда получим /2e = - • -£ • , /3e = - • -£ • .

3 ££0 3 ££0

Таким образом, локальная матрица Ke и вектор правых частей /е СЛАУ для конечного

элемента We =We + Ге имеют вышеописанный вид. Для построения (8) необходимо локальным номерам узлов 1, 2, 3 каждого элемента е поставить в соответствие глобальные номера узлов p, p = 1, P. Затем расширенные матрицы и векторы правых частей всех конечных элементов

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

Преобразуем итоговую СЛАУ с учетом краевых условий для функции ф на границе Г^ и далее

решаем ее любым из известных методов.

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

Кроме того, использование известного соотношения P =(e- 1)е0Е (Кл/м2) позволяет провести численный расчет компонент вектора поляризации в образце.

Инициализация параметров моделирования

Для запуска указанных вычислительных алгоритмов необходимо определить первоначальное распределение зарядов. Данная подзадача была решена с помощью возможностей среды визуального моделирования CASINO 2.42 с параметрами моделирования, приведенными в табл. 1.

Таблица 1

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

Название, ед. измерения Значение

Энергия падающего излучения E0, кэВ 10

Энергия останова движения Emm, кэВ 2

Количество электронов, N 1000

Данный пакет предназначен для моделирования траекторий движения электронов методом Монте-Карло [11]. Расчет проведен по отношению к модельному сегнетоэлектрическому кристаллу ниобату лития. Графическое изображение симулированного распределения траекторий электронов в образце представлено на рис. 4.

Распределение начальной плотности накопленного заряда можно описать с помощью распределения Гаусса для локальной области:

p = Poexp(- г2/52), (11)

где р0 - объемная плотность накопленного заряда, Кл/м3; г - длина радиус-вектора, проведенного в данную точку пространства из начала координат, м.

Параметры р0, 5 в (11) определяются путем прямой аппроксимации распределения, найденного с помощью метода статистических испытаний.

Рис. 4. Распределение траекторий электронов в модельном кристалле Ы№03 (стартовая энергия пучка Е0=10 кэВ, значение пороговой энергии Е =2 кэВ, N = 1000 историй электронов).

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

модели (табл. 2), а также рассчитаны физические параметры среды (табл. 3).

Таблица 2

_Физические константы, используемые в модели_

Название, ед. измерения Значение

Константа Больцмана к, Дж/К 1.38-10-23

Масса электрона т, кг 9.110-31

Заряд электрона г, Кл 1.6-10-19

Электрическая постоянная е0, Кл/(В-м) 8.85-10-12

Таблица 3

Параметры моделирования для ниобата лития_

Название, ед. измерения Значение

Коэффициент диффузии Б, м2/с 0,0275

Дрейфовая подвижность электронов тп , м2/(В-с) 1,0549

Длина свободного пробега электронов 1 , м 0.6-10-6

Температура Т, К 300

Средняя тепловая скорость движения электрона V , м/с 105

Диэлектрическая проницаемость £ (Ы№03) 30

Начальное распределение плотности р, Кл/м3 р =353ехр ~-(г2 + г2 )/(о.6 -10-6 )2"

Программная реализация модели

Программная реализация модели проведена в ППП МаАаЬ. Для проведения модельного расчета рассмотрим простейший случай - процесс релаксации заряда (в этом случае функция источника в диффузионном уравнении системы (2) отсутствует).

Результатом моделирования является поверхность плотности распределения зарядов, график которой представлен на рис. 5 (для фиксированного момента времени ^ = Ю-12 с ).

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

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

Рис. 6. Результаты моделирования напряженности поля, созданного накопленным зарядом, и профиля распределения поляризации по глубине для Ы№03 при г = 0 .

Заключение

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

1. Sogr, A.A., Maslovskaya, A.G., Kopylova, I.B. Advanced modes of imaging of ferroelectric domains in the SEM // Ferroelectrics. - 2006. - Vol. 341. - P. 29.
2. Cazaux, J. About the mechanisms of charging in EPMA, SEM, and ESEM with their time evolution // Microscopy and Microanalysis. - 2004. - Vol. 10, № 6. - P. 670-680.
3. Рау, Э.И., Евстафьева, Е.Н., Андрианов, М.В. Механизмы зарядки диэлектриков при их облучении электронными пучками средних энергий // ФТТ. - 2008. - Т. 50, № 4. - С. 599-607.
4. Кортов, В.С., Звонарев, С.В. Моделирование методом Монте-Карло транспорта электронов в заряженных при облучении кристаллических диэлектриков // Математическое моделирование. - 2009. - Т. 20, № 6. - С. 79-85.
5. Борисов, С.С., Грачев, Е.А., Негуляев, Н.Н., Черемухин, Е.А. Моделирование поляризации диэлектрика в процессе облучения электронным пучком // Прикладная физика. - 2004. - № 1 - С. 118-124.
6. Chan, DSH, Sim, KS, Phang, JCH. A simulation model for electron irradiation induced specimen charging in a scanning electron microscope // Scanning microscopy. - 1993. - Vol. 7. - P. 847-895.
7. Suga, H., Tadokoro, H., Kotera, M. A simulation of electron beam induced charging-up of insulators // Electron microscopy. - 1998. - Vol. 1. - P. 177-178.
8. Яненко, Н.Н. Метод дробных шагов решения многомерных задач математической физики -Новосибирск: Наука, 1967. - 197 c.
9. Формалев, В.Ф., Ревизников, Д. Л. Численные методы - М.: ФИЗМАТЛИТ, 2004. - 400 с.
10. Young, W.K., Bang, H. Finite Element Method Using MATLAB (Mechanical Engineering). - CRC-Press, 1996. - 544 p.
11. Drouin, D., Coutre, A.R., Joly, D., Tastet, X., Aimez, V., Gauvin, R. A fast and easy-to-use modeling tool for scanning electron microscopy and microanalysis users // Scanning. - 2007. - Vol. 29. - P. 92.
РАСТРОВОЙ ЭЛЕКТРОННОЙ МИКРОСКОПИИ
Другие работы в данной теме:
Контакты
Обратная связь
support@uchimsya.com
Учимся
Общая информация
Разделы
Тесты