Вычисление потенциала в задаче полевой эмиссии электронов из малого конического катода

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Abstract

Emission of electrons represents a phenomenon of electron emissions by solid or liquid substances under different an external, impacts: heating, sunlight, electric field etc. Field emission is one of the most significant types of emission, as potentially this emission mechanism allows the largest emission current. This paper mainly focuses on the quantum-mechanical phenomenon in a tiny dimension (nano) silicon cathode. This work shall assumingly yield a solution for equation for the potential in the case of a single emitter.



Аннотация

Электронная эмиссия - явление испускания электронов из твердых тел или жидкостей. Полевая эмиссия является одним из наиболее значимых видов эмиссии, так как потенциально этот механизм эмиссии допускает самый большой эмиссионный ток. Исследование нацелено на изучение данного квантово- механического явления в кремниевом катоде малых (нано) размеров. Ожидаемым результатом работы является решение уравнения для потенциала в случае отдельно стоящего эмиттера в вакууме.



Введение

эмиссия катод численный теплоперенос

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

Явление эмиссии известно человеку с середины 19 века. В 1839 году Александр Беккерель наблюдал фотоэффект в электролите. Позднее, в 1973 году, Уиллоуби Смит обнаружил фотопроводимость селена. Также, в 1887 году Генрих Герц установил, что при воздействии ультрафиолетом на цинковые разрядники сильно облегчается прохождение искры.

Несколько важных открытий были сделаны русским физиком Александром Столетовым в конце 19 века, в том числе был выведен первый закон внешнего фотоэффекта. Однако явление фотоэффекта было объяснено только в 1905 году Альбертом Энштейном на основе гипотезы о квантовой природе света, высказанной Максом Планком. За это открытие А. Энштейн получил нобелевскую премию в 1921 году.

Другой тип электронной эмиссии - термоэлектронная эмиссия, стал известен благодаря Фредерику Гутри в 1873 году. Он обнаружил, что нагретая до красна металлическая сфера, обладающая отрицательным зарядом, теряет его, разряжаясь в воздух. Однако при положительном заряде и аналогичных условиях разряда не происходит. В 1880 году термоэлектронная эмиссия была открыта заново Томасом Эдисоном при попытке выяснить причину разрушения нити накаливания в лампе.

Слово «Электрон» не использовалось в публикациях до 1897 года, так как он не был известен как отдельная частица. После его открытия Оуен Уилланс Ричардсон начал изучать явление, которое он в последствии и назвал «термоэлектронной эмиссией».

Явление автоэлектронной эмиссии также было открыто во второй половине 19 века. Впервые его обнаружил Р. В. Вуд в 1897 году. Попытки теоретического обоснования предпринимались В 1923 и 1928 годах Шоттки и Р. Г. Фаулером, Л. Нордгеймом соответственно.

Открытие взрывной электронной эмиссии состоялось лишь в 1966 году группой ученых из Томского института ( институт автоматизированных систем управления автоэлектронной эмиссией), Института оптики атмосферы Сибирского отделения АН СССР, а также Ленинградского государственного университета и Института автоматики и электроники.

Безусловно, теперь стоит дать определения всем известным видам электронной эмиссии: термоэлектронной, фотоэлектронной, вторичной электронной, автоэлектронной и взрывной электронной.

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

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

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

Автоэлектронная эмиссия - это явление испускания электронов проводящими электрический ток веществами под действием внешнего электрического поля без затрат энергии на предварительное возбуждение электронов вещества, что свойственно другим видам эмиссии. Это явление происходит благодаря тунелированию электронов через потенциальный барьер на поверхности вещества, что становится возможным благодаря искривлению потенциального барьера из-за воздействия сильного внешнего поля на поверхность эмиттера. Эмиссионный ток, который можно получить с помощью данного вида эмиссии, потенциально, на несколько порядков превышает токи, получаемые с помощью других видов. Другим названием этого вида эмиссии является "Холодная эмиссия". Данное явление широко используется в наноэлектронике, микроскопии, в электронной голографии и других областях.

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

В данной работе основное внимание будет уделено автоэлектронной эмиссии электронов в вакуум. Тело, из которого происходит испускание электронов, называется катодом. Для осуществления автоэлектронной эмиссии катод необходимо поместить в мощное электрическое поле с напряженностью F = 107-108 В/см. Для создания поля с такой напряженностью между макроскопическими электродами необходимо создать разность потенциалов около 10⁷вольт.

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

Стоит сказать несколько слов о механизме автоэлектронной эмиссии. Для того, чтобы электрон смог вылететь из катода, он должен преодолеть потенциальный барьер U, то есть совершить работу, направленную на преодоление сил, не позволяющих ему покинуть катод. Силы «зеркального изображения» выступают в роли основных сил, удерживающих электрон. Эти силы возникают по мере приближения электрона к поверхности катода - электрон поляризует электронный газ внутри тела, как будто создавая положительный заряд внутри катода, при чем по модулю этот заряд соответствует заряду эмиттирующего электрона. Это приводит к возникновению силы Кулона, потенциал в системе СГС имеет вид

(1.1)

где e - элементарный заряд,

x - расстояние, равное удалению свободного электрона от поверхности эмиттера.

Работа выхода, равная A = eф (ф - потенциал работы выхода), представляет собой работы, которую необходимо совершить для преодоления потенциального барьера U. Чтобы электрон преодолел потенциальный барьер, он должен обладать большей энергией чем энергия барьера. Однако, при определенных условиях, электрон способен пройти сквозь барьер не перепрыгивая его. Это становиться возможным в результате истончения потенциального барьера в следствие воздействия сильного электрического поля. Процесс прохождения электрона через потенциальный барьер называют туннельным эффектом. Именно этот эффект лежит в основе холодной эмиссии.

Суммарный потенциал в таком случае имеет вид:

(1.2)

Как было сказано выше, форма потенциального барьера изменяется - по мере усиления электрического поля, потенциальный барьер истончается. Для описания поведения потенциального барьера вводится характеристика D - прозрачность барьера. Прозрачность барьера, это безразмерная величина, находящаяся в пределах от 0 до 1, ее физический смысл - вероятность прохождения электрона сквозь потенциальный барьер в вакуум, при попадании его на границу катода. Прозрачность потенциального барьера, согласно квантово-механическим расчетам, может быть вычислена по формуле (1.3):

(1.3)

где h - постоянная Планка,

m - масса частицы,

U - потенциальная энергия,

E - энергия туннелирующей частицы.

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



1. Актуальность темы

В настоящий момент, явление автоэлектронной эмиссии слабо изучено, однако это направление является перспективным. Данное явление лежит в основе технологии создания дисплеев. Они носят названия FED (Field Emission Display) или SED (Surface conduction electron Emitter display)[11]. Данная технология позволяет получать плоские экраны с большой диагональю и рядом преимуществ по сравнению с распространенными ныне LED-дисплеями. Принцип работы FED-дисплея несколько напоминает работу экрана с кинескопом. В последнем, электронно-лучевая трубка создает поток электронов, падающий на поверхность, покрытую люминофором - экран. Отклоняемый с помощью магнитного поля поток проходит построчно по всем точкам экрана, тем самым создавая кадр. Однако, для эффективного управления лучом электронов требуется значительное пространство, вследствие чего ЭЛТ-телевизоры имеют большие габариты и вес. Технология FED-дисплеев лишена этого недостатка. Дисплей состоит из пары прозрачных панелей, разнесенных на несколько миллиметров друг от друга. На одной поверхности находятся излучатели электронов, на другой люминофор, аналогичный используемому в ЭЛТ. В отличие от кинескопа, за каждый пиксель отвечает свой излучатель, а точнее, отдельный для каждого из трех цветов - красного, синего, зеленого. Как можно было заметить, из-за отсутствия необходимости пробегать по всему экрану лучом кинескопа, размеры конструкции существенно меньше. Наряду с этим, из-за отсутствия развертки изображения - изображение создается целыми кадрами за раз, потенциальная частота обновления изображения составляет порядка 240 Гц, что больше, чем у представителей лучших ЖК-дисплеев. Еще одним преимуществом данной технологии является «живучесть» дисплеев - существуют образцы, которые при выходе из строя до 20% излучателей не имеют на дисплее «мертвых» пикселей. Вместе с тем, FED-дисплеи, наряду с ЭЛТ-дисплеями обладают сочным и красочным изображением, отменным отображением черного цвета и контрастностью до 100 000:1, что значительно превышает показатели LED-дисплеев. Время отклика изображения при использовании данной технологии достигает 1 мс, а углы обзора стремятся к теоретическому максимуму - 180 градусов, так как экран сам является источником света, а не проецирует свет сквозь себя.

Рис. 1.Принцип работы FED-панели

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

Одной из наиболее значимых работ в области изучения автоэлектронной эмиссии, является работа Р. Г. Фаулера и Л. Нордгейма. Если известно количество электронов, попадающих на барьер (что может быть рассчитано из теории твердого тела), то становится возможным рассчитать эмиссионный ток j электронов, проходящих в вакуум, и получить формулу для тока автоэлектронной эмиссии. Необходимые квантово-механические расчеты впервые были выполнены в работе [1]:

(2.1)



где y = 3,79*10^-4(F^1/2)/ф, входящие функции табулированы.

Функция t(y), находящаяся в предэкспоненциальном множителе, близка к единице и слабо изменяется с изменением аргумента.

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

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

При помещении катода в сильное электрическое поле, в нем происходят различные тепловые процессы. Помимо тривиальных явлений, стоит упомянуть некоторые неординарные. Таковым является эффект Ноттингама - выделение теплоты на катоде при автоэлектронной эмиссии и ее поглощение при термоавтоэлектронной эмиссии, обусловленные разницей между средней энергией электронов, приближающихся к поверхности электрода и покидающих его. При низких температурах, то есть автоэлектронной эмиссии, распределение электронов по энергиям близко к распределению Ферми при абсолютном нуле и через потенциальный барьер U проходят электроны с энергией несколько меньшей уровня Ферми, при этом катод нагревается за счет энергии электронов, приходящих на свободные уровни из электрической цепи. При термоэлектронной эмиссии (в случае высокой температуры) электроны покидают эмиттер с уровней выше, чем уровень Ферми. Заполнение этих уровней электронами из цепи приводит к охлаждению катода. Данное явление открыто У. Б. Ноттингамом в 1941 году [2].

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

В статье Данилова В. Г., Руднева В. Ю., Гайдукова Р. К, Кретова В. И. [4] рассматривается метод моделирования теплопереноса в автоэмиссионном нанокатоде. В основе лежит модифицированная задача Стефана с условиями Ноттингама на вершине эмиттера и его свободной границе. В численном моделировании используется модифицированная система фазового поля. С помощью численного моделирования исследуется влияние явления Ноттингама на динамику распространения границы раздела жидкой и твердой фаз в эмиттере.

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

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

Изображение прототипа модели представлено на рис. 2.1, геометрические размеры в таблице 1. Модель, использованная при расчетах представлена на рисунке 2.2



Таблица 1

Высота катода	10-15 мкм
Диаметр основания катода	6 мкм
Радиус скругления вершины катода	15 нм
Угол при вершине катода	20о

Рис. 1.1

Рис. 1.2



Пространственные переменные меняются в следующих пределах:

На рисунке обозначены границы r₁(t) и r₂(t), которые являются свободными и возникают в результате проплавления катода - между r₁(t) и r₂(t) катод находится в жидком состоянии, через знак ? обозначена боковая поверхность. Область, в которой происходит моделирование, представляет собой усеченный конус, где основания задаются как

r = R₀ и r = R.

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

В сферических координатах оператор Лапласа имеет вид

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



Таблица 2

Разогрев катода при прохождении через него электрического тока описывается уравнением (2.2) - уравнением теплопроводности.

(2.2)

Температура обезразмеривается по формуле (2.3)

(2.3)

Из следующих выражений (2.4) определяются безразмерное время и радиус (координата)

(2.4)

Из предыдущего получается уравнение теплопроводности в сферических координатах:

(2.5)

Удельная проводимость определяется формулой (2.6), плотность же тока j_in внутри катода находится из уравнения (2.7) - уравнения неразрывности.

(2.6)

(2.7)

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

Отдельно стоит сказать про условие на верхнем основании катода (острие), оно определяется эффектом Ноттингама, которое является условием охлаждения. Краевое условие, определяющее краевой баланс на острие катода, имеет вид (2.8)

(2.8)



Где j_em - плотность эмиссионного тока, ?_Nott - средняя энергия эмитирующих электронов, Т - реальная температура по Кельвину, ? - реальная координата в метрах.

В безразмерной форме краевое условие примет вид (2.9)

(2.9)

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

Уравнение для потенциала получается из (2.7) и имеет вид:

(2.10)

Где краевыми условиями являются:

Таблица 3

На нижней границе	U0
На боковых гранях
На верхней границе	U1 = 0

2. Постановка задачи

Настоящая задача несколько отличается от рассмотренной ранее. В данной задаче изменения касаются самой модели, еще более приближенной к настоящим условиям. Геометрически, новая модель представляет собой усеченный конус, на верхнем основании которого (острие) находится полусфера, радиус которой равен радиусу верхнего основания конуса. При решении задачи рассматривается сечение, принадлежащее плоскости, проходящее через ось описанного тела вращения. Катод находится в прямоугольной области, при чем нижнее основание катода совпадает с границей области, а верхняя граница, являющаяся анодом, отстоит на некоторое малое расстояние от катода. Модель так же имеет две границы r₁(t) и r₂(t), которые являются свободными и возникают в результате проплавления катода. На рисунке 3.1 представлена данная модель (рис. 3.2 демонстрирует вершину катода в той же модели).

Рис. 3.1

Рис. 3.2

За основу взяты размеры из таблицы 1. Таблица 4 содержит дополнительную информацию.



Таблица 4

Высота катода	10.015 мкм
Диаметр основания катода	6 мкм
Радиус скругления вершины катода	15 нм
Высота границы r1	h1
Высота границы r2	h2
Высота области	10.03 мкм
Ширина области	6.03 мкм

В ходе моделирования h₁ и h₂ принимают различные значения.

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

Уравнение для потенциала имеет вид:

(3.1)

V - объем всей области, где решается задача. Будем полагать, что коэффициент является константой во всей расчетной области. Умножим (3.1) на пробную функцию , и проинтегрируем по частям:

(3.2)

На границе области полагаем = 0, поэтому (3.2) примет вид

(3.3)

Если (3.3) верно для любой функции , то тогда (3.3) принимают за определение обобщенного решения уравнения (3.1), в том числе и в случае разрывного коэффициента .

В рассматриваемой работе прямоугольная область V разбивается на две подобласти: V_ext - объем снаружи катода, V_int - объем внутри катода c разными значениями разрывного коэффициента Интеграл в левой части (3.3) в этом случае разбивается на два:

(3.4)

Подставим (3.4) в (3.3), применим тождество (3.5) и формулу Остроградского-Гаусса (3.6) - математическую формулу, которая выражает поток непрерывно-дифференцируемого векторного поля (А) через замкнутую поверхность интегралом от дивергенции этого поля по объёму, ограниченному этой поверхностью:

(3.5)

(3.6)

(3.7)

в V_ext и на границе:

в V_int и на границе:



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

(3.8)

Положим и

, (3.9)

где Г_S - боковая поверхность катода, L - верхняя полусфера (поверхность эмиссии). Тогда из (3.3) остались интегралы

, (3.10)

На поверхности (кривой) L справедливо равенство

, (3.11)

где D - средний коэффициент прозрачности барьера. Тогда



(3.12)

Катод окружен вакуумом, поэтому можно считать равным 0. Тогда

(3.13)

- дельта-функция.

Краевое условие на верхней границе катода для уравнения теплопроводности

, (3.14)

- напряженность поля. k - коэффициент, полученный путем объединения коэффициентов из уравнения (2.9). На верхней границе . Если распределение температуры известно, то можно получить условие вида

. (3.15)

Краевые условия для уравнения для потенциала приведены ниже:



Таблица 5

На нижней границе	U0
На боковых границах
На верхней границе	U1 = 0

3. Методы

Для решения поставленной задачи необходимо выбрать один из нескольких методов решения задач данного класса:

1. Метод конечных разностей (МКР)

2. Метод конечных элементов (МКЭ)

Метод конечных разностей -- численный метод решения дифференциальных уравнений, основанный на замене производных разностными схемами. Является сеточным методом.

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

Далее, не вдаваясь глубоко в описание методов, будет приведено лишь их сравнение:

Таблица 6

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



Таким образом, МКЭ является наиболее оптимальным для решения данной задачи потому, что:

· Область, в которой решается задача является геометрически сложной.

· В задаче используется разрывный коэффициент

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

Метод конечных элементов

Реализация метода состоит из нескольких шагов:

1. Триангуляция

Пусть область Щ ограничена, - кусочно-гладкая замкнутая кривая. Разобьем область на треугольники , i = 1,2,...,n. Такое разбиение называется триангуляцией. Треугольники, у которых две вершины находятся на границе области, могут иметь криволинейную сторону. Эту сторону заменим отрезком прямой и получившуюся область обозначим через . Вершины треугольников называются узлами триангуляции, занумеруем их и пусть D1 ? множество узлов, принадлежащих границе, D2 ?множество внутренних узлов,

12 D = D1 D2 ? множество всех узлов триангуляции, а J1, J2, J ? количество узлов в D1, D2, D соответственно.

2. Построение конечных элементов

Рассмотрим любой треугольник ? . Для каждой вершины этого треугольника определим функцию , равную единице в этой вершине и равную нулю в других вершинах (здесь i - номер треугольника, j - номер вершины в множестве всех узлов D). Простейшей такой функцией является линейная функция вида z = ax + by + c, где z = 1 в одной из вершин треугольника и z = 0 в двух других вершинах.

Пример. Пусть, например, треугольник Д7 имеет вершины P2(-1;2), P5(1;0), P11(3;4). Тогда для него определены три функции , , верхний индекс обозначает номер треугольника, а нижний - номер вершины треугольника (узла триангуляции). Здесь функция задает плоскость, проходящую через точки (1;2;1),(1;0;0),(3;4;0), функция задает плоскость, проходящую через точки (1;2;0),(1;0;1),(3;4;0), а функция задает плоскость, проходящую через точки (1;2;0),(1;0;0),(3;4;1). Выпишем эти функции, используя стандартное уравнение плоскости, проходящей через три данные точки. После простых преобразований получаем: Итак, каждому треугольнику Д_i соответствуют три линейные функции , , где i- номер треугольника, а j1, j2, j3 - номера его вершин.

Сопоставим j-му узлу произвольный параметр , j = 1, …, J. Пусть вершины треугольника Д_i имеют в множестве D номера j1, j2, j3. Определяем новую функцию

(4.1)

В нашем примере получаем

где коэффициенты a, b, c вычисляются по формулам:



Из построения функций следует, что значение в узле с номером j равно . Функция называется конечным элементом.

3. Построение системы уравнений

Определим в функцию . Из определения видно, что внутри каждого треугольника функция дифференцируема, так как совпадает с линейной функцией . Кроме того, . Если треугольники Д_i1 и Д_i2 имеют общую сторону или вершину, то на этой стороне или вершине . Поэтому непрерывна и кусочно-дифференцируема на всей области Если узел с номером j принадлежит границе области, то, учитывая граничное условие, положим

(4.2)

Тогда функция равна 0 на границе области и. Ее графиком является поверхность, составленная из плоских треугольников. В определении функции участвуют J переменных . Из них J₁ переменных определяются граничным условием (4.2). Требуется выбрать значения остальных J₂ переменных так, чтобы функция была в каком-то смысле близка к решению исходной краевой задачи. Рассмотрим выражение

(4.3)

Очевидно, I_n(u_n) - функция от J₂ переменных . Индекс n в I_n указывает на отличие от функционала I, так как в (4.3) область интегрирования , а не .

Из необходимого условия экстремума получаем для всех внутренних узлов

(4.4)

Система (4.4) линейна, так как I_n(u_n) - квадратичная функция. Порядок системы равен J₂. Ввиду очевидной связи между функцией I_n(u_n) и функционалом I(u) естественно предположить, что

4. Структура матрицы системы (4.4).

Запишем систему (4.4) в виде

(4.5)

где

(4.6)

Заметим, что присутствует только в тех слагаемых суммы (4.6), для которых область интегрирования, то есть Д_i, содержит j - й узел. Поэтому после дифференцирования по в сумме (4.6) исчезнут все слагаемые, не содержащие . Значит, каждое уравнение системы (4.4) содержит только относительно небольшую часть от общего числа неизвестных. Поэтому большая часть элементов матрицы B равна 0. Такие матрицы называются разреженными (sparse). Для решения систем с разреженными матрицами разработаны специальные методы.

5. Вычисление элементов матрицы системы (4.4).

Рассмотрим одно слагаемое в (4.6). Имеем

Учитывая линейность функций , запишем их в виде

Отсюда получаем: и

, .

(4.7)

(4.8)



Кроме того, из (4.1) следует, что Заметим, что

Поэтому с учетом формулы (4.8) получим

(4.9)

Аналогично вычисляются производные по от всех слагаемых в сумме (4.6), в которых используется узел с номером j₁. То же рассуждение верно и для любого узла. Элементы b_j,k матрицы B из уравнения (4.5) и f_j вектора F вычисляем следующим образом:

· j - е уравнение системы получаем, дифференцируя функционал I_n(u_n) по переменной и используя формулы (4.6) и (4.9);

· в полученном уравнении находим слагаемые, содержащие ;

· элемент b_j,k равен сумме коэффициентов при ;

элемент fj равен сумме всех интегралов вида входящих в j- е уравнение.

После того, как матрица B и вектор F вычислены, решение системы осуществляется с помощью стандартных компьютерных процедур.

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

Выбор пал на MATLAB (сокращение от англ. «Matrix Laboratory», в русском языке произносится как Матламб) -- пакет прикладных программ для решения задач технических вычислений и одноимённый язык программирования, используемый в этом пакете. Конкретнее на пакет PDETool Box.

Таким образом, исходную задачу можно разделить на три этапа:

1. Триангуляция области

2. Вывод системы линейных алгебраических уравнений (СЛАУ)

3. Решение уравнения

Триангуляция области будет осуществляться отдельно, 2 и 3 пункт программируются в MATLAB.

4. Расчетная область

Подготовка расчетной области происходила в несколько этапов:

1. Построение прямоугольной сетки

2. Построение триангуляции

3. Уточнение геометрии области

Построение прямоугольной сетки (преимущественно прямоугольной) реализовано с помощью языка программирования Python 3.0. Была написан следующий код:

import matplotlib.pyplot as plt

import numpy as np

def setKa(H,D,r,step,step2):

# H-высота катода в мкм

# D-диаметр катода в мкм

# r-радиусскругления вершины катода в нм

# step-основной шаг сетки

# step2-дополнительный шаг сетки(шаг на полукруге)

X = np.array([-500*D,-r,r,500*D])

Y = np.array([0,1000*H,1000*H,0])

F1 = np.arange(0, np.pi*(2*r+1)/(2*r), np.pi/(2*r))

F2 = np.arange(0, np.pi, np.pi/(2*r))

fig, ax = plt.subplots()

ax.plot(X, Y, color = "blue")

ax.plot(r*np.cos(F1),1000*H + r*np.sin(F1), color = "blue")

xx = np.arange(-625*D,(625+2*step)*D,step)

yy = np.arange(0,1100*H+step,step)

i = 0

while i <= 7500/step:

if ((xx[i] <= -r)|(xx[i] >= r)):

ax.plot([xx[i],xx[i]], [0,1100*H], color = "green", lw=0.2)

else:

ax.plot([xx[i],xx[i]], [0,1000*H], color = "green", lw=0.2)

ax.plot([xx[i],xx[i]], [1000*H+r,1100*H], color = "green", lw=0.2)

i+=1

i = 0

while i <= 11000/step:

if ((yy[i] <= 1000*H)|(yy[i] >= 1000*H+r)):

ax.plot([-625*D,625*D],[yy[i],yy[i]], color="green",lw=0.2)

else:

ax.plot([-625*D,-r], [yy[i],yy[i]], color = "green", lw=0.2)

ax.plot([r,625*D], [yy[i],yy[i]], color = "green", lw=0.2)

i+=1

i = 0

while i < r:

i+=step

ax.plot(i*np.cos(F1),1000*H + i*np.sin(F1), color = "green", lw=0.2)

i = 0

while i <= 2*r:

ax.plot([r*np.cos(F1[i]),0],[1000*H + r*np.sin(F1[i]),1000*H], color = "green", lw=0.2)

i+=step2

plt.xlim((-2*r,2*r))

plt.ylim((1000*H-2*r,1000*H+2*r))

#plt.xlim((-4000,4000))

#plt.ylim((0,11000))

ax.set_aspect('equal')

plt.show()

результатом выполнения кода является сетка, изображенная на рисунке 5.1 (часть сетки).

Рис. 5.1

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

Вторая версия устранила имевшиеся недостатки за счет неравномерности прямоугольной сетки и корректировки условий для окрестности полукруга. Пример кода ниже, результат на рисунке 5.2:

import matplotlib.pyplot as plt

import numpy as np

def setKa(H,D,r,step,step2,step3):

# H-высота катода в мкм

# D-диаметр катода в мкм

# r-радиусскругления вершины катода в нм

# step-основной шаг сетки

# step2-дополнительный шаг сетки(шаг на полукруге)

X = np.array([-500*D,-r,r,500*D])

Y = np.array([0,1000*H,1000*H,0])

F1 = np.arange(0, np.pi*(2*r+1)/(2*r), np.pi/(2*r))

F2 = np.arange(0, np.pi, np.pi/(2*r))

fig, ax = plt.subplots()

ax.plot(X, Y, color = "blue")

ax.plot(r*np.cos(F1),1000*H + r*np.sin(F1), color = "blue")

xx1 = np.arange(-625*D,(625+step)*D,step)

xx2 = np.arange(-625*D,(625+step3)*D,step3)

yy1 = np.arange(0,1000*H-r+step,step)

yy2 = np.arange(1000*H-r,1000*H+2*r+step3,step3)

i = 0

while i <= 7500/step:

ax.plot([xx1[i],xx1[i]], [0,1000*H-r], color = "green", lw=0.2)

i+=1

i = 0

while i <= 7500/step3:

if ((xx2[i] <= -r)|(xx2[i] >= r)):

ax.plot([xx2[i],xx2[i]], [1000*H-r,1000*H+2*r], color = "green", lw=0.2)

else:

ax.plot([xx2[i],xx2[i]], [1000*H-r,1000*H], color = "green", lw=0.2)

ax.plot([xx2[i],xx2[i]], [1000*H+r*(1-(xx2[i]/r)**2)**(1/2),1000*H+2*r], color = "green", lw=0.2)

i+=1

i = 0

while i <= (1000*H-r)/step:

ax.plot([-625*D,625*D],[yy1[i],yy1[i]], color="green",lw=0.2)

i+=1

#i = 0

i = 0

while i <= (3*r)/step3:

if ((yy2[i] <= 1000*H)|(yy2[i] >= 1000*H+r)):

ax.plot([-625*D,625*D],[yy2[i],yy2[i]], color="green",lw=0.2)

else:

ax.plot([-625*D,-r*(1-((yy2[i]-1000*H)/r)**2)**(1/2)], [yy2[i],yy2[i]], color = "green", lw=0.2)

ax.plot([r*(1-((yy2[i]-1000*H)/r)**2)**(1/2),625*D], [yy2[i],yy2[i]], color = "green", lw=0.2)

i+=1

i = 0

while i < (r-step2):

i+=step2

ax.plot(i*np.cos(F1),1000*H + i*np.sin(F1), color = "green", lw=0.2)

i = 0

while i <= 2*r:

ax.plot([r*np.cos(F1[i]),0],[1000*H + r*np.sin(F1[i]),1000*H], color = "green", lw=0.2)

i+=step2

plt.xlim((-2*r,2*r))

plt.ylim((1000*H-2*r,1000*H+2*r))

#plt.xlim((-4000,4000))

#plt.ylim((0,11000))

ax.set_aspect('equal')

plt.show()

Рис. 5.2

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

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

Геометрия области задается с помощью матрицы g, имеющей размерность 12хМ, где М - число сегментов области. Общая структура имеет вид: каждый столбец определяет 1 сегмент. Первый элемент в столбце определяет тип сегмента (2 -- отрезок прямой, 4 -- эллипса); 2, 3 (4, 5) строки определяют x (y) координаты начальной и конечной точки соответственно; 6, 7 строки -- номер области слева и справа по направлению обхода. Для отрезков прямых следующие строки не нужны и задаются нулями; для эллипса 8, 9 строки определяют x и y координаты центра эллипса, а 10, 11 строки -- его x и y полуоси (для окружности они совпадают и равны радиусу окружности); 12 строка определяет угол поворота эллипса вокруг центра против часовой стрелки (в радианах).

Изначально, в Матлаб была перенесена расчетная область, использовавшаяся при создании прямоугольной сетки. Матрица описания области таблица 8:

Таблица 8

2	2	2	2	2	2	2	4	4	2
0,3	-0,3	-0,3015	-0,3015	0,3015	0,3015	0,3	0,0015	0	-0,0015
-0,3	-0,3015	-0,3015	0,3015	0,3015	0,3	0,0015	0	-0,0015	-0,3
0	0	0	1,003	1,003	0	0	1	1,0015	1
0	0	1,003	1,003	0	0	1	1,0015	1	0
0	0	0	0	0	0	2	2	2	2
2	1	1	1	1	1	1	1	1	1
0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	1	1	0
0	0	0	0	0	0	0	0,0015	0,0015	0
0	0	0	0	0	0	0	0,0015	0,0015	0
0	0	0	0	0	0	0	0	0	0

Внешний вид области (цифрами обозначены соответствующие номера сегментов) рис. 5.3, 5.4 (увеличенное изображение).

Рис. 5.3

Рис. 5.4

Триангуляция осуществлялась с посредством встроенной функции, ее синтаксис:

[p, e, t] = initmesh( g, `Hmax', 0.0015),

Где параметр g - описание геометрии области, `Hmax' - установка ограничения размера треугольника при разбиении, “0.0015” - максимальный размер треугольника.

Выходные параметры, p - массив узлов конечноэлементной сетки (столбцам соответствуют узлы):

- первая строка - горизонтальные координаты узлов

- вторая строка -вертикальные координаты узлов;

e - матрица граничных элементов на границах раздела зон:

столбцам соответствуют граничные элементы (стороны конечных элементов, принадлежащие границам раздела зон или границе расчётной области);

первые две строки - номера начальных и конечных узлов граничных элементов;

строки 3, 4 - длина «дуги» от начала граничного сегмента до начального и конечного узла граничного элемента, отнесённая к длине «дуги» граничного сегмента;

строка 5 - номера граничных сегментов, которым принадлежат граничные элементы;

строки 6, 7 - номера зон, примыкающих слева и справа к граничным элементам;

t - матрица треугольных конечных элементов (столбцам соответствуют треугольники):

- t(1:3,ie) - глобальные номера узлов треугольника с номером ie,

- t(4,ie) - номер зоны, которой принадлежит треугольник с номером ie.



Рис. 5.5

Внешний вид триангуляции (части) c размером элемента, сопоставимым радиусу полуокружности, представлен на рисунке 5.5.

Данная сетка включает 464715 треугольников. На ее построение Вычислительная машина оснащенная четырехъядерным процесором по 3.5 Ггц каждое и двухпоточной памятью 2х4 Гб затратила более 5 часов. Полученная сетка слишком крупная для провидения вычислений с достаточной точностью, а получение сетки с размером треугольника хотя бы в 3 раза меньшим данного займет более суток, не говоря уже о провидении дальнейших вычислений. Было принято решение предоставить возможность алгоритму пакета построить сетку в адаптивном режиме c большим размером треугольника (0.15), рисунок 5.6.

Рис. 5.6

Очевидно, что данная сетка непригодна для проведения вычислений по той же причине - крупность треугольника, но видна особенность - в местах сближения сегментов области сетка уточняется самостоятельно.

Пакет PDETool содержит функцию, позволяющую переопределить сетку с уточнением на основе существующей, ее вид:

[p1, e1, t1] = refinemesh(g, p, e, t),

Где p1, e1, t1 - новые массивы данных после уточнения, соответствующие p, e, t.

Уточненный вариант области представлен на рисунке 5.7.

Рис. 5.7

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

Геометрия области изменена следующим образом:

· Добавлена “расплавленная” часть катода;

· добавлена вертикальная полоса симметричная относительно прямой Х=0, шириной 2r (r - радиус полукруга), в связи с особенностями правой части уравнения.

Эти изменения привели к увеличению матрицы g с 12х10 до 32х10, и число внутренних подобластей увеличилось с 2 до 12. Ниже представлены обновленная матрица g (таблица 9) и ее визуализация (рисунки 5.8, 5.9 - увеличенная часть катода). Используя данную матрицу, была получена триангуляция области, которая после однократного уточнения приняла вид, представленный на рисунках 5.10,5.11 (увеличенная часть катода).

Таблица 9

2	2	2	2	2	2	2	2	2	4	4	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2
0,0015	-0,3	-0,3015	-0,0015	0,3015	0,3015	0,3	0,00254475	0,00194775	0,0015	0	-0,0015	-0,00194775	-0,00254475	0,0015	-0,0015	-0,0015	0,3	-0,3015	0,0015	0,00194775	-0,0015	-0,00254475	0,0015	0,0015	-0,0015	0,0015	-0,0015	0,0015	-0,0015	0,0015	-0,0015
-0,0015	-0,3015	-0,3015	0,0015	0,3015	0,3	0,00254475	0,00194775	0,0015	0	-0,0015	-0,00194775	-0,00254475	-0,3	-0,0015	0,0015	-0,3	0,0015	-0,0015	0,3015	0,0015	-0,00194775	-0,0015	0,00254475	0,0015	-0,0015	0,0015	-0,0015	0,0015	-0,0015	0,0015	-0,0015
0	0	0	1,003	1,003	0	0	0,9965	0,9985	1	1,0015	1	0,9985	0,9965	0,9985	0,9965	0	0	1,003	1,003	0,9985	0,9985	0,9965	0,9965	0	0	0,9965	0,9965	0,9985	0,9985	1	1
0	0	1,003	1,003	0	0	0,9965	0,9985	1	1,0015	1	0,9985	0,9965	0	0,9985	0,9965	0	0	1,003	1,003	0,9985	0,9985	0,9965	0,9965	0,9965	0,9965	0,9985	0,9985	1	1	1,003	1,003
0	0	0	0	0	0	3	5	7	11	11	11	6	4	10	10	0	0	0	0	5	6	6	5	8	4	10	6	11	8	12	2
9	2	2	12	1	1	1	1	1	12	12	2	2	2	11	9	4	3	2	1	7	8	4	3	3	8	5	10	7	11	1	12
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0,0015	0,0015	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0,0015	0,0015	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0

Рис. 5.8

Рис. 5.9



Рис. 5.10

Рис. 5.11

5. Результаты вычислений

Используя полученное ранее описание геометрии области, триангуляцию, а также вспомогательные матрицы, полученные вместе с триангуляцией области, с помощью Matlab было получено распределение функции (напряжения) внутри расчетной области, в том числе и внутри катода.

Уравнение вида - div(c*grad(U)) + a*U = f может быть решено с помощью функции U = assempde( b, p, e, t, c, a, f ) пакета PDETool Matlab. Где U - распределение искомой величины в центрах треугольниках сетки , b - матрица краевых условий для всех сегментов, согласованная с матрицей описания геометрии g, она имеет следующую структуру:

· 1 строка - внешний или внутренний сегмент;

· 2 строка - условие Дирихле или Неймана;

· 3 строка - длина выражения для q;

· 4 строка - длина выражения для g;

· 5 строка - длина выражения для h;

· 6 строка - длина выражения для r;

· 7 и далее строки - ASCII-коды q, g, h, r,

где q, g, h, r это коэффициенты в одном из видов условий: либо h*U = r, либо c*grad(U) + q*U = g; p, e, t - массивы полученные при триангуляции ранее; c - строка, содержащая значения коэффициента «» для каждой подобласти, а - строка, содержащая значения коэффициента «а» для каждой подобласти, f - строка, содержащая значения правой части для каждой подобласти.

Правая часть уравнения имеет вид

- уравнение полуокружности, верхней границы эмиттера, [a, b] = [-0.0015, 0.0015], экспонента здесь выступает в роли дельта-функции.

При r₁ = 0.9985, r₂ = 0.9965 , a = 0 во всей области, с = 0.001 снаружи катода, 1000 внутри катода, 10000 в расплавленном слое, получены следующие результаты (рисунки 6.1, 6.2, 6.3):

Рис. 6.1 (плоскость ХУ)

Рис. 6.2 (плоскость ХУ)



Рис. 6.3 (трехмерное изображение с сеткой)

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

При r₁ = 0.9985, r₂ = 0.9965 , a = 0 во всей области, с = 0.001 снаружи катода, 1000 внутри катода, 10000 в расплавленном слое, получены следующие результаты (6.4, 6.5, 6.6):



Рис. 6.4 (плоскость ХУ)

Рис. 6.5 (трехмерное изображение)



Рис. 6.6 (трехмерное изображение с сеткой)

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

Напряженность можно вычислить с помощью встроенной функции

.

Графическое представление напряженности по Uy в первом и во втором случаях представлены на рисунках 6.7 и 6.8 соответственно.



Рис. 6.7

Рис. 6.7

Напряженность по OY постоянна внутри катода, в вакууме по мере увеличения координаты Y уменьшается с небольшой скоростью в большей части расчетной области на обоих рисунках, изменения наблюдаются при приближении к верхней границе расчетной области - внутри катода напряженность постоянна, вне катода резко уменьшается по мере приближения к "острию" катода, при чем в первом случае падение происходит до нескольких сотен, а во втором до нескольких тысяч (со знаком минус). И в том, и в другом случае максимальная по модулю напряженность располагается на оси симметрии катода (х = 0), выше катода.



Список литературы

1. Fowler R. H., Nordheim L. Electron emission in intense electric fields // Proceedings of the Royal Society of London. Series A. - 1928. - Vol.119, no. 781. P. 173 -181

2. Nottingham W. В. Remarks on Energy Losses Attending Thermionic Emission of Electrons from Metals // Physical Review. - 1941.

3. Furcey G. Field Emission in Vacuum Microelectronics. - Springer, 2005. - 205 p.

4. Данилов В. Г., Руднев В. Ю., Гайдуков Р. К., Кретов В. И. Математическое описание режима «плавление затвердевание» в автоэмиссионном катоде при учете эффекта Ноттингама // Наноструктуры. Математическая физика и моделирование - 2013. Т. 9. - С. 39-84.

5. Данилов В. Г., Руднев В. Ю., Гайдуков Р. К., Кретов В. И. Математическое моделирование эмиссии из катодов малых размеров // 2014. 232с.

6. Дюжев Н. А., Гудкова С. А., А. Махиборода М., Федирко В. А. Исследование эмиссионных свойств кремниевых катодов различной геометрии // Вакуумная наука и техника, материал 12 научно-технической конференции с участием зарубежных специалистов / Под ред. Д. В. Быкова. - М.: МИЭМ, 2005. - С. 221-224.

7. Г. Стренг, Дж. Фикс. Теория метода конечных элементов. // МИР - 1977

8. Matlab. Partial Differential Equation Toolbox™ User's Guide // 2019

9. Доля П.Г. Использование MATLAB. Решение дифференциальных уравнений...