Содержание к диссертации
Введение
Глава 1. Комбинированный метод конечных элементов- интегрального уравнения для расчета свойств диэлектрических волноводов вблизи критических частот 19
1.1. Постановка задачи 20
1.1.1. Внешняя краевая задача 24
1.1.2. Численные методы для расчета направляемых мод ДВ в приближении слабой волноводности 29
1.2. Основы метода 42
1.2.1. Метод конечных элементов 44
1.2.2. Метод интегрального уравнения 49
1.2.3. Объединенное матричное уравнение 53
1.3.1 Ірограммная реализация метода 56
1 .3.1. Задание геометрии системы 56
1.3.2. Разбиение области определения решения на элементы (триангуляция) 57
1.3.3. Формирований матриц и решение проблемы собственных значений 59
1.3.4. Вывод и оценка результатов 61
1.4. Тестовые расчеты 61
1.4.1. Круглые и эллиптические волноводы 64
1.4.2. Волноводы прямоугольного сечения 68
1.5. Расчет диэлектрических волноводов с вытянутым поперечным сечением 71
1.5.1. Дисперсионные зависимости ДВ с вытянутым поперечным сечением 76
1.5.2. Критические частоты ДВ с вытянутым поперечным сечением 79
1.6. Расчет волноводов с поперечным сечением невыпуклой формы 83
1.7. Выводы к главе 1 87
Глава 2. Исследование вытекающих мод диэлектрических волноводов сложного поперечного сечения 89
2.1. Постановка задачи 90
2.1.1. Свойства вытекающих мод ДВ 90
2.1.2. Идентификация направляемых и вытекающих мод ДВ с потерями или активной средой 94
2.1.3. Внешняя краевая задача 101
2.1.4. Обзор существующих численных методов 102
2.2. Модификация метода для случая вытекающих мод 106
2.2.1. Модификация уравнений метода 107
2.2.2. Определение вида функции Грина 110
2.2.3. Итерационная процедура 1 13
2.2.4. Особенности программной реализации метода 1 15
2.3. Результаты тестовых расчетов 117
2.3.1. Вытекающие моды круглого ДВ 119
2.3.2. Вытекающие моды эллиптического ДВ 121
2.3.3. Расчет вытекающих мод ДВ типа «канал» 123
2.3.4. Расчет вытекающих мод эллиптического ДВ типа «канал» 125
2.4. Исследование высекающих мод многослойного ДВ типа «трубка» 127
2.4.1. Многослойный ДВ типа «трубка» (брэгговское волокно) 130
2.4.2. Влияние дефекта типа «щель» на свойства основной вытекающей моды многослойного ДВ типа «трубка» 1 33
2.4.3. Влияние дефекта в виде неконцентричиости 137
2.4.4. Влияние эллиптической деформации многослойного ДВ типа «трубка» на свойства основной вытекающей моды 140
2.5. Выводы к главе 2 146
Глава 3. Развитие методов расчета устройств вакуумной микроэлек троники, содержащих конструктивные элементы с сильно различаю щимися масштабами 148
3.1. Постановка задачи 149
3.1.1. Внутренняя краевая задача 152
3.1.2. Обзор существующих численных методов для расчета МАЭК 154
3.2. Основы метода 160
3.2.1. Метод конечных элементов 162
3.2.2. Метод интегрального уравнения 165
3.2.3. Объединенное матричное уравнение 174
3.3. Результаты тестовых расчетов 176
3.4. Расчет матричного автоэмиссионного катода 180
3.5. Расчет эмиссионных свойств пористого кремния 1 84
3.6. Моделирование автоэлектронной эмиссии с многомасштабной (фрактальной) поверхности 191
3.6.1. Множество Жюлиа как модель фрактальной поверхности 192
3.6.2. Методика численного моделирования автоэлсктронной эмиссии с фрактальной поверхности 195
3.6.3. Расчет вольтамперных характеристик 198
3.6.4. Влияние фрактальной размерности на эмиссионные характеристики 198
3.7. Выводы к главе 3 203
Заключение 204
Литература 208
- Численные методы для расчета направляемых мод ДВ в приближении слабой волноводности
- Разбиение области определения решения на элементы (триангуляция)
- Расчет волноводов с поперечным сечением невыпуклой формы
- Особенности программной реализации метода
Введение к работе
Одним из наиболее интенсивно развивающихся в последние два десятилетия разделов радиофизики является вычислительная электродинамика -наука, занимающаяся разработкой алгоритмов и методик моделирования электромагнитных полей [1-9].
Совершенствование элементной базы современных ЭВМ, а также разработка новых методов численного решения уравнений в частных произвол-пых, позволяет в настоящее время решать широкий класс электродинамических задач, таких, как расчет СВЧ резонаторов и волноведуших структур, устройств интегральной оптики, замедляющих систем для СВЧ усилителей и генераторов, элементов электронных ускорителей, антенных систем, и так далее, без каких—либо дополнительных предположений и упрощений, исходя непосредственно из уравнений Максвелла и материальных уравнений, описывающих свойства заполняющих систему сред.
Среди большого числа численных методов, используемых для этого, следует особенно выделить 3 метода, обладающих особой мощностью и универсальностью и, по этой причине, наиболее популярных. Это метод конечных разностей (в последнее время особенно бурно развивается вариант это і-о метода во временном представлении - Finite Difference Time Domain method [1-3]), метод конечных элементов (МКЭ) [4,5] и метод интегрального уравнения (МИУ) [7,8]. Последний метод в электродинамике также часто называют методом граничных элементов.
Тем не менее, в электродинамике существует значительное число задач, которые не поддаются адекватному рассмотрению с помощью стандартных вариантов упомянутых методик. Это обстоятельство связано, главным образом, с появлением все новых устройств и систем, как в СВЧ, так и в оптическом диапазоне, которые, во-первых, могут обладать достаточно сложной геометрической формой, и, во-вторых, работать в новых, ранее не используемых режимах и областях управляющих параметров.
В частности, во многих СВЧ устройствах, представляющих собой неоднородные в поперечном сечении линии передачи (микрополосковые антенны, диэлектрические волноводы (ДВ) различного типа, щелевые линии и т.д. [10-13]) в качестве рабочих мод используются направляемые волны (в частности, в ДВ - поверхностные моды), которые формально могут распространяться без затухания, если не учитывать потери в среде. Однако в ряде случаев эти системы работают в диапазоне частот, непосредственно примыкающем к критической частоте направляемой моды, при этом размер области, занимаемой полем моды в поперечном направлении, значительно возрастает и может превышать поперечные размеры самой сердцевины волновода во много раз. Точно на частоте отсечки можно считать, что поперечный размер ноля становится бесконечным. Отметим, что задача определения критических частот направляемых мод ДВ сама по себе является чрезвычайно важной, так как этим определяется, например, область одпомодового режима работы волновода. Для решения этой задачи также необходимо учитывать бесконечный, по сути, размер поля в поперечном сечении, то есть то обстоятельство, что ДВ является открытой электродинамической структурой. Ниже критической частоты направляемая мода ДВ перестает существовать и распространение волны вдоль системы сопровождается излучением ее энергии через боковые сгенки волновода. Таким образом, можно выделить две особенности рассматриваемых электродинамических систем, которые существенно затрудняют их теоретический анализ и практический расчет:
Эти системы имеют сложную геометрическую форму в поперечном сечении и сложный закон распределения диэлектрической проницаемости по поперечному сечению.
Эти системы являются "открытыми" электродинамическими структурами, что приводит к необходимости учитывать значительную протяженность нолей собственных и несобственных мод в поперечном сечении, а также возможное излучение энергии из волновода через его боковую поверхность.
Похожие проблемы возникают и при исследовании таких новых оптических систем, как фотонные кристаллы [14-16], полые и микроструктурированные волокна [17-20], представляющих собой кварцевую или стеклянную микроструктуру с периодическими повторяющимися воздушными капиллярами, многослойные (брэгговские) волоконные световоды.
Наличие периодичности в направлении, поперечном к направлению распространения волны, приводит к тому, что такая структура представляет собой, фактически, многомерную периодическую систему [21]. Дефект структуры, соответствующий отсутствию одного или нескольких воздушных отверстий, играет роль, аналогичную роли сердцевины в обычных волоконных световодах. В этом случае возникает область с повышенным значением эффективного показателя преломления, и в результате интерференции падающих и отраженных от такой неоднородности волн формируется направляемая мода микроструктурированного волокна.
Другим механизмом, обеспечивающим канализацию энергии вдоль волокон рассматриваемого типа, является наличие фотонных запрещеЕшых зон (photonic gap), или, на языке высокочастотной электродинамики, зон непропускания [22]. Фотонная запрещенная зона обеспечивает высоких коэффициент отражения для излучения, распространяющегося вдоль полой сердцевины волокна, что позволяет существенно уменьшить оптические потери, присущие модам полых волноводов [19,20]. По существу, на том же самом принципе работают многослойные (брэгговские) волоконные световоды [23,24].
Важным обстоятельством, отличающим как микроструктурированные
волокна с полой сердцевиной, так и брэгговские волокна, является то, что они всегда обладают конечными (хотя и малыми) потерями на излучение, даже если считать, что материал, из которого изготовлен световод, имеет нулевой тангенс потерь. Это излучение объясняется просачиванием энергии водны через запрещенные области за счет туннельного эффекта, так как в реальных системах периодическая структура всегда конечна в поперечном на-
правлении. Тем не менее, за счет выбора подходящих значений параметров, эти потери могут быть сделаны чрезвычайно малыми. Следовательно, эти системы являются принципиально «открытыми», то есть их электродинамические характеристики определяются частичным излучением электромагнитных волн через боковую поверхность волновода при одновременной канализации основной энергии волны вдоль его оси.
Таким образом, актуальной является задача расчета собственных и несобственных (вытекающих) мод открытых электродинамических волноводов (диэлектрических волноводов и оптических волоконных световодов), имеющих сложную форму поперечного сечения с неоднородным в этом сечении распределением диэлектрической проницаемости.
Основными методами, используемыми в данной диссертационной работе для решения этой задачи, являются метод конечных элементов и метод интегрального уравнения . Для обоих этих методов характерна высокая универсальность, то есть возможность решать широкий класс задач с помощью единого подхода. Оба метода позволяют использовать относительно гладкую аппроксимацию границ расчетных областей и учитывать сложное распределение диэлектрической проницаемости.
Стандартный вариант МКЭ [4,5,25,26] описывает распределение поля только внутри ограниченной области. Разработаны его модификации, позволяющие учитывать поля в открытом пространстве [4], например, применяя бесконечные элементы или импедансные граничные условия. Существенным ограничением таких подходов является использование аппроксимаций, описывающих затухание поля с удалением от волновода, с помощью параметров, которые необходимо подбирать отдельно. К сильным сторонам МКЭ можно отнести простоту и эффективность описания сложной геометрии системы и свойств материалов.
1 В работе; не рассматривается метол конечных разностей, ідавньїм образом, из-за существенных технических трудностей, вочникаюших и ) гоч методе при ею компьютерной реализации, применительно к расчету систем со сложной криволинейной формой границы.
В МИУ для описания поля в открытой части пространства используется только его значения на некотором контуре или в некоторой конечной области. Таким образом, МИУ позволяет рассчитывать широкий класс открытых волноводов вне зависимости от поперечного размера поля. Однако, описание сложного распределения диэлектрической проницаемости и диэлектрических потерь в данном случае оказывается сложной задачей.
Проведенный анализ показывает, что объединение методов МКЭ и МИУ для расчета открытых структур может позволить использовать их сильные стороны и обойти недостатки. Комбинированный метод конечных элементов - интегрального уравнения способен эффективно описывать как распределение полей в открытом пространстве, так и сложную структуру сердцевины волновода или волоконного световода.
Следует отметить, что комбинированный метод конечных элементов -интегрального уравнения используется в электродинамике уже достаточно давно [см. например, 4, гл. 9 и 27], однако объектом подобных исследований являлись, главным образом, задачи рассеяния электромагнитных волн на.не-однородностях сложной формы. Применительно к задачам расчета собственных направляемых мод открытых диэлектрических волноводов комбинированный метод развивался только в работе [28], однако развиваемый в ней вариант МКЭ не учитывает трудности, возникающие в конечно-разностных методиках при моделировании векторных полей, в частности, возможность появления паразитных решений [29]. Кроме того, используемые в работе [28] алгоритмы позволяют практически исследовать только системы с достаточно простой геометрией.
Что касается использования комбинированного метода для расчета несобственных (в частности, вытекающих) волн открытых диэлектрических волноводов СВЧ и оптического диапазона, то такие работы до настоящего времени отсутствовали.
Комбинированный метод КЭ-ИУ оказался плодотворным подходом также при решении другой актуальной проблемы, возникающей в современ-
ной физической электронике — разработке методов моделирования электрических полей в автоэмиссионпых устройствах вакуумной микроэлектроники [30,31]. При расчете автоэмиссионных катодов возникают следующие сложности:
Большой диапазон масштабов. Если расстояние катод-анод в таких системах имеет порядок 10 Зм, то радиус скруглення острия эмиттера может составлять порядка 10~н м. При этом неоднородности и микроструктура поверхности на острие эмиттера, сильно влияющие на эмиссионные свойства системы, еще на 1-2 порядка меньше.
Максимальная точность расчета необходима именно на микроостриях, имеющих минимальный размер.
В отличие от задачи расчета открытых волноводов, в данном случае МИУ используется не для описания полей в открытой части пространства, а для высокоточного описания поля на эмитирующих остриях.
Несмотря на существенные отличия в характере задачи и необходимость использовать несколько отличные модификации МКЭ и МИУ, идеи и алгоритмы сопряжения используемых методов близки как для задачи расчета собственных мод открытых волноводов, так и при расчете устройств вакуумной микроэлектроники. Это обстоятельство позволяет в рамках данной работы получать решения обеих перечисленных задач с использованием единой численной методики.
Целью настоящей диссертации является исследование направляемых и вытекающих мод открытых диэлектрических волноводов со сложной формой поперечного сечения и сложным распределением диэлектрической проницаемости в поперечном сечении а также автоэмиссионных систем вакуумной микроэлектроники, имеющих многомасштабную структуру.
Для достижения поставленной цели в диссертации решаются следующие задачи.
Разрабатывается комбинированный метод конечных элементов -интегрального уравнения для расчета направляемых мод открытых диэлектрических волноводов с неоднородной сердцевиной сложной формы, в том числе возле частоты отсечки и непосредственно на критической частоте.
Комбинированный метод обобщается для расчета вытекающих мод открытых диэлектрических волноводов с неоднородной сердцевиной сложной формы, с учетом потерь или усиления в среде.
На основе предложенных алгоритмов разрабатывается и тестируется комплекс программ для ЭВМ, предназначенных для расчета электродинамических параметров (дисперсии, затухания за счет излучения, распределения полей собственных и несобственных мод) открытых диэлектрических волноводов и оптических световодов.
С использованием созданных программ проводятся расчеты дисперсионных характеристик и полей собственных поверхностных мод диэлектрических волноводов со сложной формой поперечного сечения. Особое внимание уделяется исследованию указанных величин вблизи критических частот и точно на критических частотах направляемых мод.
Проводится исследование диэлектрических волноводов с вытянутой формой поперечного сечения (с большим отношением поперечных размеров в двух перпендикулярных направлениях).
Изучается влияние разнообразных дефектов на дисперсию и потери практически значимых систем оптического диапазона (трехслойных брэгговских волокон), работающих на вытекающих модах.
Развивается вариант комбинированного метода конечных элементов-интегрального уравнения для расчета автоэмиссионных устройств вакуумной микроэлектроники.
Проводится проверка пределов применимости приближенной аналитической модели, используемой для оперативных расчетов параметров ячеек матричных автоэмиссионных катодов.
Решается задача расчета напряженности поля на поверхности эмиссионных центров пористого кремния, поверхность которого представляет собой упорядоченную структуру периодически расположенных высоких острий.
Проводятся исследования многомасштабной микроструктуры поверхности эмиттера, аппроксимируемой фрактальным множеством Жюлиа.
Научная новизна результатов работы состоит в следующем.
Развит новый универсальный численный метод для расчета направляемых мод, в том числе на критической частоте, а также вытекающих мод широкого класса открытых диэлектрических волноводов и волоконных световодов, имеющих сложную форму поперечного сечения, неоднородный профиль диэлектрической проницаемости в сердцевине волновода и, возможно, потери или усиление в среде.
Впервые проведено исследование трансформации электродинамических параметров диэлектрических волноводов и многослойных брэг-говских волокон с вытянутой формой поперечного сечения при постепенном увеличении параметра формата.
.3. Впервые ггроведено исследование влияния геометрических дефектов трехслойных брэгговских оптических волокон на дисперсионные характеристики и параметры затухания вытекающих мод. Результаты расчетов показывают, что дефекты типа "щель" и "неконцентричность" оказывают незначительное влияние на свойства волокон (в частности, на дисперсию и потери на излучение) в широком .диапазоне изменения параметров дефектов. Эллиптическая деформация трехслойного брэг-
говского волокна оказывает незначительное влияние на дисперсию, но приводит к значительному уменьшению потерь на излучение. 4. Комбинированный метод конечных элементов - интегрального уравнения впервые применен для расчета устройств вакуумной микроэлектроники, что позволило обеспечить высокую точность расчета напряженности электрического поля на поверхности эмитирующего острия.
Практическая значимость работы. Развитые в диссертационной работе методы расчета открытых диэлектрических волноводов могут быть использованы при разработке новых СВЧ - устройств (диэлектрических антенн, сенсорных датчиков, направленных ответвителей на связанных диэлектрических волноводах). Разработанные методики применимы также для расчета ряда перспективных устройств волоконной оптики, таких как фотонные кристаллы, микроструктурированные и брэгговские волоконные световоды. Предложенные методики допускают обобщение на более сложный класс открытых электродинамических систем, в частности, диэлекрических волноводов с большим относительным изменением показателя преломления по сечению, волноводов с оболочкой в виде плоско-слоистой среды и т.д.
Результаты исследования эффектов усиления электрического поля на эмитирующих поверхностях в присутствии микронеоднородностей представляют интерес для создания перспективных устройств вакуумной микроэлектроники, использующих новые материалы для автоэмиссионных катодов (пористый кремний, углеродные пленки).
В диссертации разработаны следующие программные комплексы для решения практически важных задач при помощи комбинированного метода: 1. Пакет прикладных программ для расчета направляемых и вытекающих мод диэлектрических волноводов и оптических световодов в приближении слабонаправляемых мод. Программы обладают значительной универсальностью по отношению к геометрической форме исследуе-
мых систем, распределению относительной диэлектрической проницаемости и диэлектрических потерь в сердцевине, позволяют учитывать бесконечный размер полей в открытом пространстве и обеспечивают расчет параметров волноводов во всем диапазоне частот. 2. Пакет прикладных программ для расчета распределения полей в устройствах вакуумной микроэлектроники. Программы обладают универсальностью по отношению к геометрической форме исследуемых систем, распределению относительной диэлектрической проницаемости в пределах ячейки автоэмиссиошюго катода, позволяют добиться значительной точности в определении напряженности поля на поверхности эмиттера.
Программы использовались при выполнении исследований по г/б НИР 'ГОР, проводимых в 2001-2004 гг. в Саратовском государственном университете, грантов РФФИ №02-02-17317 и № 03 -02-16161.
Основные положения и результаты, выносимые на защиту.
Комбинированный метод конечных элементов - интегрального уравнения позволяет проводить эффективный расчет собственных (направляемых) и несобственных (вытекающих) мод открытых электродинамических систем - диэлектрических волноводов и оптических световодов - во всем диапазоне частот, включая малую окрестность критической частоты, вблизи которой поля мод неограниченны в пространстве.
Для диэлектрических волноводов с сильно вытянутой формой поперечного сечения существуют две области частот с различным характерным поведением направляемых собственных волн. Для частот, больших некоторого граничного значения, дисперсионные характеристики и пространственные распределения полей мод близки к характеристикам и распределениям полей, свойственным двумерным (плоским) диэлектрическим структурам. Для частот, меньших граничного
значения, дисперсия и структура полей вытянутых волноводов соответствуют этим характеристикам в трехмерных системах. Граничное значение частоты зависит от параметра формата волновода (соотношения поперечных размеров по двум перпендикулярным направлениям).
Геометрические дефекты трехслойных брэгговских волокон в виде продольной щели или неконцентрического смещения слоев оказывают незначительное влияние на дисперсию и потери на излучение вытекающих мод из волокон, в широком диапазоне изменения параметров дефектов. В то же время, эллиптическая деформация волокна в поперечном сечении приводит к заметному снижению потерь на излучение.
Для автоэмиссионного катода, состоящего из упорядоченной структуры периодически расположенных высоких острий, при фиксированном расстоянии между ними, существует оптимальное расстояние катод-анод, для которого коэффициент усиления поля на поверхности эмиттера достигает максимального значения.
Личный вклад соискателя.
В работах [160-162,165-167,169-172] автору принадлежит вывод соотношений, разработка алгоритмов и программ, реализующих комбинированный метод КЭ - ИУ, проведение расчетов, участие в анализе и интерпретации результатов.
В работах [156-159,163,164,168] вклад автора заключается в разработке алгоритмов и программ расчета потенциала, создаваемого проводником, в форме множества Жюлиа, напряженности поля и плотности тока эмиссии вдоль заданной эквипотенциали, а также в проведении расчетов при помощи этих программ.
Апробация работы и публикации.
Материалы диссертационной работы докладывались на 5-й Международной школе по хаотическим колебаниям CHAOS'98 (Саратов, 1998г), на школе-семинаре «Нелинейные дни в Саратове для молодых» (Саратов, 1998г), на ХІ-й и ХП-й зимних школах по СВЧ электронике и радиофизике (Саратов, 1999 и 2003 г.г.), на Международной конференции «Электроника и радиофизика СВЧ UHF'99» (С.-Петербург, 1999 г.), на 12-й Международной конференции по вакуумной микроэлектронике (Дармштадт, Германия, 1999 г.), на 5-м и 7-м рабочих семинарах ШЕЕ Penza-Saratov Chapter (Саратов, 2000 и 2002 г.г.), на Международной межвузовской конференции «Современные проблемы электроники и радиофизики СВЧ (Саратов, 2001 г.), на 7-й Всероссийской конференции студентов-физиков ВНКСФ-7 (С.Петербург, 2001 г.), на 4-й Международной конференции ТЕБЕ по вакуумным источникам электронов (TVECS-2002), (Саратов, 2002 г.), на 6-м Международном симпозиуме по электрическим и магнитным полям (EMF-2003) (Аахен, Германия, 2003 г.).
Материалы диссертации обсуждались на научных семинарах кафедр электроники, колебаний и волн и нелинейной физики СГУ. По теме диссертации опубликовано 3 статьи в научных журналах, 6 статей в трудах научных конференций и 8 - тезисов докладов.
Структура и объем работы
Диссертация состоит из введения, трех глав, заключения и списка цитированной литературы. Диссертация содержит 226 страниц текста, 67 рисунков и графиков и список литературы из 172 наименований.
Краткое содержание работы
Первая глава посвящена решению актуальных задач радиофизики, касающихся расчета свойств направляемых мод диэлектрических волноводов, в
том числе вблизи и на частоте отсечки. Проведено исследование трансформации дисперсионных характеристик поверхностных мод открытых ДВ с вытянутым поперечным сечением вблизи критической частоты при постепенном увеличении отношения поперечных размеров волновода в двух перпендикулярных направлениях.
Показано, что для ДВ с сильно вытянутой формой поперечного сечения существуют две области частот с различным характерным поведением направляемых собственных волн. Для частот, больших некоторого граничного значения, дисперсионные характеристики и пространственные распределения полей мод близки к характеристикам и распределениям полей, свойственным двумерным (плоским) диэлектрическим структурам. Для частот, меньших граничного значения, дисперсия и структура полей вытянутых волноводов соответствуют этим характеристикам в трехмерных системах.
Получена и подтверждена численно оценка граничного значения безразмерной частоты, разделяющей области различного поведения направляемых мод. Граничное значение зависит от отношения сторон сечения волновода. Данная оценка позволяет определить, можно ли использовать для описания конкретного ДВ модель планарного волновода, либо необходимо использовать трехмерную модель ДВ.
Для высших направляемых мод с сильно вытянутой формой поперечного сечения исследована применимость планарного приближения непосредственно на частоте отсечки. Рассчитаны критические частоты ДВ с поперечным сечением невыпуклой формы (в форме двулистника). Подобная геометрия хорошо моделирует поперечное сечение направленного ответвителя, широко используемого в интегральной оптике.
Для расчета свойств направляемых мод диэлектрических волноводов развит комбинированный метод конечных элементов - интегрального уравнения. Метод обеспечивает расчет ДВ с произвольной формой сердцевины и произвольным распределением диэлектрической проницаемости по сечению волновода.
Во второй главе рассматриваются вопросы расчета вытекающих мод диэлектрических волноводов, для чего обобщается комбинированный метод конечных элементов - интегрального уравнения, развитый в первой главе диссертации. Проведены исследования практически важных волоконно-оптических устройств, работающих на вытекающих модах трехслойных брэгговских волокон.
Впервые проведено исследование влияния геометрических дефектов трехслойных брэгговских оптических волокон на дисперсионные характеристики и параметры затухания вытекающих мод. Показано, что дефекты "щель" и "неконцентричность" оказывают незначительное влияние на свойства волокон (в частности, на дисперсию и потери на излучение) в широком диапазоне изменения параметров дефектов. Эллиптическая деформация трехслойного брзгговского волокна оказывает слабое влияние на дисперсию, но приводит к заметному уменьшению потерь на излучение. Дастся качественное объяснение этому эффекту.
Третья глава посвяшена решению актуальных задач физической электроники - расчету и моделированию ряда автоэмиссионных устройств, применяемых в вакуумной микроэлектронике. Определена область применимости приближенной аналитической модели ячейки матричного автоэмиссионного катода, предложенной для оперативных расчетов напряженности поля на острие эмиттера, в зависимости от параметров ячеек. Рассчитана напряженность поля на поверхности эмиссионных центров пористого кремния, представляющих собой плотно упакованную периодическую структуру близко расположенных высоких острий. Показано, что данной системы возможно так подобрать геометрические параметры, чтобы достичь максимального значения коэффициента усиления поля на поверхности. Проведено исследование автоэлектронной эмиссии с самоподобной (фрактальной) поверхности.
Для моделирования устройств вакуумной микроэлектроники разработана модификация комбинированного метода конечных элементов - интегрального уравнения.
Численные методы для расчета направляемых мод ДВ в приближении слабой волноводности
В электродинамике разработано большое количество численных методов, среди которых наиболее широко известны метод конечных элементов (МКЭ), метод конечных разностей и метод интегрального уравнения (МИУ). Однако краевая задача (1.23)-(1.26) обладает рядом специфических особенностей, которые должны обязательно учитываться при использовании той или иной методики. Во-первых, вблизи от частоты отсечки поперечный размер поля направляемой моды резко возрастает и на частоте отсечки стремится к бесконечности. Таким образом, метод должен учитывать распределение полей в открытом пространстве. Во-вторых, в последние годы кроме круглых диэлектрических волноводов с постоянным значением относительной диэлектрической проницаемости є широко исследуются многослойные диэлектрические волноводы (ДВ), ДВ с прямоугольным, эллиптическим, и более сложными формами сечения, а также различные волокна с переменным распределением диэлектрической проницаемости є по сечению. В этом разделе приведен обзор существующих численных методов, применимых к решению краевой задачи (1.23)-(1.26). Особый акцент сделан на возможность применения этих методов для расчета ДВ сложной формы поперечного сечения, переменным распределением к по сечению и в области частот около отсечки. В конце обзора будет приведена сравнительная оценка сильных и слабых сторон каждого метода. Метод интегрального уравнения. Начнем с метода интегрального уравнения по площади. Метод интегрального уравнения по площади является универсальным методом и может применяться для расчета ДВ с переменным профилем є, и даже для расчета анизотропных систем (в векторном варианте). Будем рассматривать простейшую геометрию, когда относительная диэлектрическая проницаемость є(х,у) является непрерывной функцией поперечных координат внутри конечной односвязной области Qr, граница которой - контур Л . Удобным для приложения является интегральное уравнение [37], которое в скалярном приближении выглядит так [38]: где -относительная диэлектрическая проницаемость вне конту Макдональда. Алгебраизация уравнения (1.27) может быть произведена различными способами. Удачным является метод коллокаций, который можно рассматривать как вариант метода Галеркина [39]. При таком подходе выбирается определенная система узлов [29] с координатами г: и на решение налагается условие обращения в нуль невязки в этих узлах. Здесь f - система базисных функций. В подробной записи эти соотношения имеют вид однородной системы линейных алгебраических уравнений с матрицей ДДи ) отн области частот около отсечки. В конце обзора будет приведена сравнительная оценка сильных и слабых сторон каждого метода. Метод интегрального уравнения.
Начнем с метода интегрального уравнения по площади. Метод интегрального уравнения по площади является универсальным методом и может применяться для расчета ДВ с переменным профилем є, и даже для расчета анизотропных систем (в векторном варианте). Будем рассматривать простейшую геометрию, когда относительная диэлектрическая проницаемость є(х,у) является непрерывной функцией поперечных координат внутри конечной односвязной области Qr, граница которой - контур Л . Удобным для приложения является интегральное уравнение [37], которое в скалярном приближении выглядит так [38]: где -относительная диэлектрическая проницаемость вне конту Макдональда. Алгебраизация уравнения (1.27) может быть произведена различными способами. Удачным является метод коллокаций, который можно рассматривать как вариант метода Галеркина [39]. При таком подходе выбирается определенная система узлов [29] с координатами г: и на решение налагается условие обращения в нуль невязки в этих узлах. Здесь f - система базисных функций. В подробной записи эти соотношения имеют вид однородной системы линейных алгебраических уравнений с матрицей ДДи ) относительно неизвестных коэффициентов С . Из условий существования нетривиального решения этой системы следует характеристическое (нелинейное) уравнение для поперечного волнового числа w. Недостаток этого подхода в том, что существуют определенные ограничения на выбор узлов, при нарушении которых вычислительный процесс может расходиться. Анализ этих ограничений представляет собой дополнительную задачу, которую во многих случаях удается решить [37,38,40]. Согласно результатам [29], для круглых диэлектрических волноводов метод интегрального уравнения по площади обеспечивает практически ту же высокую точность, что и метод обыкновенных дифференциальных уравнений, являющийся оптимальным для расчета круглых волокон [41]. Однако, метод интегрального уравнения по площади почти на порядок медленнее. Отметим, что метод интегрального уравнения по площади позволяет быстро рассчитать поле в волокне при использовании разложения по полиномам Чебышева [29]. Коэффициенты разложения легко определяются методом обратных осительно неизвестных коэффициентов С . Из условий существования нетривиального решения этой системы следует характеристическое (нелинейное) уравнение для поперечного волнового числа w. Недостаток этого подхода в том, что существуют определенные ограничения на выбор узлов, при нарушении которых вычислительный процесс может расходиться. Анализ этих ограничений представляет собой дополнительную задачу, которую во многих случаях удается решить [37,38,40]. Согласно результатам [29], для круглых диэлектрических волноводов метод интегрального уравнения по площади обеспечивает практически ту же высокую точность, что и метод обыкновенных дифференциальных уравнений, являющийся оптимальным для расчета круглых волокон [41]. Однако, метод интегрального уравнения по площади почти на порядок медленнее. Отметим, что метод интегрального уравнения по площади позволяет быстро рассчитать поле в волокне при использовании разложения по полиномам Чебышева [29]. Коэффициенты разложения легко определяются методом обратных итераций [42]. Заметим также, что скалярную задачу можно свести к обычной задаче на собственные значения (когда собственное значение входит в уравнение линейно), если считать заданным числом р, а неизвестным числом к. В этом случае время расчета сокращается в несколько раз, однако для практических приложений такой подход менее удобен. Для расчета ДВ с гладкой границей Л можно использовать систему метагармонических функций [43] где J№ -функции Бесселя. Как уже отмечено выше, при использовании метода коллокаций существуют некоторые ограничения на выбор узлов. В работе [37]
Разбиение области определения решения на элементы (триангуляция)
Назначение препроцессора - активное взаимодействие с пользователем. Основной его целью является обеспечение комфортной работы пользователя и реализация необходимого числа сервисных функций, для того, чтобы обеспечить решение задач различных видов при различных граничных условиях. Препроцессор должен иметь стандартный интерфейс и достаточный набор инструментов для корректного описания возможных вариантов входных данных и методов решения задачи. Важным является и то, что программа должна быть удобна в освоении и использовании и доступна для пользователя, не обладающего специальными знаниями в области численных методов. Была разработана программа-препроцессор, удовлетворяющая вышеперечисленным требованиям. Препроцессор имеет стандартный пользовательский интерфейс и написан на языке Delphi. Разработанная программа предоставляет следующие возможности: - позволяет задать геометрию системы, определить области с различным значением относительной диэлектрической проницаемости є, - задавать свойства генерируемой сетки - расположение областей сгущения и разрежения сетки, разбиение границы сопряжения Г, - сохранять и редактировать данные в режиме, удобном для пользователя, просматривать статистику, работать одновременно с произвольным количеством открытых файлов и т.п. - предоставляет дополнительные инструменты - генерация часто встречающихся геометрических фигур. В верхней части рисунка 1.8 изображен общий вид интерфейса пользователя разработанной программы. 1.3.2. Рачбиение области определения решения на элементы (триангуляция) В ходе работы препроцессора пользователь вводит необходимую информацию о свойствах сетки. Необходимо разбить область определения решения на треугольные элементы, учитывая разнообразные требования, заданные пользователем, например: - максимальная площадь треугольного элемента; - области сгущения и разрежения сетки; - требования к качеству треугольных элементов (минимальный внутренний угол) и др. - области с различными значениями диэлектрической проницаемости. Нужно заметить, что от качества сгенерированной сетки существенно зависит точность полученного решения. Для генерации сетки используется мощная программа Triangle [78], полученная из общедоступной библиотеки программного обеспечения NETLIB [79]. Программа Triangle написана на стандартном С и может быть откомпилирована на любой системе, где существует стандартный компилятор С. На входе Triangle получает файл с информацией о системе, созданный препроцессором.
На выходе программа выдает файлы с данными триангуляции. Triangle также может использоваться для улучшения уже сгенерированной сетки. В ходе работы системы программа Triangle вызывается из препроцессора, что облегчает работу с этим мощным и сложным генератором сетки. В нижней части рисунка 1.8 показан результат работы программы Triangle. Препроцессор используется для просмотра сгенерированной сетки и управления процессом триангуляции (изменение параметров, итеративное сгущение сетки и так далее). 1.3.3. Формирований матриц и решение проблемы собственных значений После того, как входные данные, включая конечноэлементную сетку, подготовлены и записаны в файл, следующей стадией решения является формирований матриц и репгение результирующего матричного уравнения (1.53). Для этого была разработана программа-решатель СОММА. Программа написана на языке Фортран с использованием нескольких внешних библиотек. Для заданного пользователем параметра w, матричные элементы матриц A(w2) и В из (1.53) формируются согласно формулам, приведенным выше (раздел 1.2). Матрицы хранятся в сжатом формате хранения разреженных матриц, описанном в [80]. Для решения обобщенной алгебраической проблемы собственных значений матричного уравнения (1.53) используется метод Арнольди [81], реализованный в пакете программ ARPACK[82]. Метод Арнольди эффективен для нахождения собственных чисел и собственных векторов несимметричных разреженных матриц, его суть состоит в вычислении ортогонального базиса Qk подпространства Крылова, являющегося линейной оболочкой векторов {q,Cq,Cq2 ,...,СкAq], q- произвольный начальный вектор, С = В1 A(w) (см.уравнение(1.53)). Собственные значения матрицы Нк Q[CQk порядка к хороню приближают крайние числа спектра матрицы С уже при небольших к . Для решения частной задачи решения системы нелинейных уравнений с квадратной несимметричной разреженной матрицей используется специализированный пакет UMF [83]. Программа базируется на многофронтальном методе факторизации больших разреженных линейных несимметричных матриц [84]. В результате решения обобщенной алгебраической проблемы собственных значений мы получаем спектр собственных чисел к и соответствующие собственные вектора, состоящие из найденных значений напряженности поля в узлах сетки и нормальной производной поля в узлах сетки, лежащих на границе сопряжения. Таким образом, для заданного внешнего поперечного числа w рассчитаны частоты собственных мод (так как к = о?іс) и распределения напряженности поля для этих мод. Задавая на входе программы w — 0, можно рассчитать частоты отсечки для заданных собственных мод. На практике следует вводить в программу число, близкое к нулю (например, w- 10 so) так как при w= 0 функция Мак-дональда К{) в функции Грина (1.43) обращается в бесконечность. В программе реализована также возможность решения задачи расчета поперечного волнового числа и для заданной частоты. Так как эта процедура необходима, в частности, для расчета вытекающих мод ДВ, то она подробно описана в следующей главе. Отметим достоинства предлагаемой в данной диссертации методики в сравнении с работой [28], где также используется комбинация МКЭ и интегрального уравнения. Отличие настоящей работы от [28] состоит в применении более мощных и надежных алгоритмов, позволяющих добиваться большей точности. Например, в работе [28] алгоритм дискретизации задачи требует однородной прямоугольной сетки (при этом прямоугольник делится на два треугольных элемента). В настоящей работе благодаря использованию мощной профаммы -триангулятора возможно произвольное сгущение сетки в области больших градиентов поля с целью увеличения точности. Также
Расчет волноводов с поперечным сечением невыпуклой формы
В цилиндрической системе координат (г,ср) граница сердцевины описывается формулой где p arctg(y/x), г2 = х1 + у1. В этой формуле параметры а и а определяют размер и форму кривой. При a - 0 двулистник переходит в окружность радиуса а. При a — 1 (a 1) такая структура может служить моделью двух соприкасающихся оптических волокон. Подобная геометрия хорошо моделирует поперечное сечение направленного ответвителя, широко используемого в интегральной оптике. На рис. 1.17 приведены сечения рассматриваемого ДВ при двух промежуточных значениях а, а именно от = 0.1 и а = 0.7. На этом рисунке изображены также многоугольники, которые окружают сердцевину ДВ и являются границами сопряжения Г. В таблице 1.4 приведены рассчитанные значения критической частоты Vc первых высших мод ДВ описанной формы при нескольких значениях параметра а. Обозначения мод были те же, что и выше. Безразмерную частоту определяли как V - kl у]г - ге где 1у — максимальная ордината двулистника в используемой системе координат (см. рис. 1.17). Как и выше, предполагалось, что диэлектрическая проницаемости сердцевины ДВ постоянна и равна є. Интересно отметить, что при увеличении параметра а изменяется порядок следования мод, если их упорядочить по возрастанию соответствующих значений Ус. В частности, при малых а второй и третьей высшими молами являются волны LPl2 и LP , тогда как при ог = 0.7 они меняются местами, т.е. происходит пересечение кривых Vc(a). На рис. 1.18 представлены линии уровня модуля поля Е(х,_у) моды /.Я, для ДВ с сечением в виде двулистника на частоте отсечки (Vt. =1.9287). Рисунок построен для случая, когда параметр а = 0.5. Заметим, что при а « ] структура мод ДВ с рассматриваемым сечением во многом похожа на структуру мод двух близко расположенных связанных волноводов[34, 55]. Классификация мод, введенная выше, в таких системах становится, как правило, в значительной мере условной. Например, распределение поля вдоль оси Ох для моды LPn при а — 0.5 имеет вид двугорбной кривой, т.е. формально первый индекс для этой моды то должен равняться трем (см.рис.1.19). Для нее было сохранено указанное выше обозначение с т-1, так как при а 0 она переходит в моду с одним экстремумом поля вдоль этой оси. 1. Выполнен критический анализ универсальных численных методов, применяемых для расчета полей направляемых мод открытых диэлектриче ских волноводов в скалярном приближении.
Обзор показал, что не существу ет достаточно универсального и удобного метода, который позволял бы с высокой точностью рассчитывать параметры направляемых мод диэлектрических волноводов со сложной формой поперечного сечения и произвольным распределением диэлектрической проницаемости в сердцевине, во всем диапазоне частот, где мода существует, в том числе и на частоте отсечки 2. Разработан комбинированный метод конечных элементов - инте грального уравнения для расчета полей направляемых мод ДВ в скалярном приближении. Метод обеспечивает расчет ДВ с произвольной формой серд цевины и произвольным распределением диэлектрической проницаемости по сечению волновода. Метод пригоден для расчетов вблизи и на частоте отсеч ки. Предлагаемая методика может быть обобщена для исследования диэлек трических волноводов более сложной конструкции: состоящих из материалов с комплексной диэлектрической проницаемостью (потери или активная сре да), имеющих оболочку в виде плоско-слоистой среды и т.д. 3. Разработан программный комплекс, реализующий предложенный метод расчета диэлектрических волноводов. Сопоставление результатов тестовых расчетов с данными, полученными другими методиками, показали точность, надежность, высокую эффективность и универсальность метода. 4. Проведены расчеты дисперсии направляемых мод и критических частот для практически важных конструкций диэлектрических волноводов со сложной формой сердцевины (эллиптической, прямоугольный ДВ, волновод в форме двулистника). 5. Исследована трансформация дисперсии направляемых мод и их собственных полей при изменении параметра формата в диапазоне / = 1 -100. 6. Для диэлектрических волноводов большого формата, работающих вблизи критической частоты, получена и подтверждена оценка граничного значения частоты, выше которой ДВ ведет себя подобно двумерной системе, а ниже - как трехмерная структура. Граничное значение зависит от формата волновода. Данная оценка позволяет определить, можно ли использовать для описания конкретного ДВ модель планарного волновода, либо необходимо использовать трехмерную модель диэлектрического волновода. Глава 2. Исследование вытекающих мод диэлектрических волноводов сложного поперечного сечения Настоящая глава посвящена вопросам разработки комбинированного метода конечных элементов - интегрального уравнения для расчета полей вытекающих мод диэлектрического волновода (ДВ). В начале главы представляется используемая модель ДВ, описаны свойства вытекающих мод, формулируется внешняя краевая задача. Далее производится обзор существующих численных методов, применимых к решению данной задачи. В результате произведенного анализа делается вывод о целесообразности приложения комбинированного метода конечных элементов - интегрального уравнения для расчета вытекающим мод ДВ. Далее приведены основные уравнения метода, описаны и обоснованы изменения, внесенные в метод, особенности формирования и решения матричных уравнений, программной реализации метода. Развитый метод тестируется при расчетах вытекающим мод круглого ДВ, для которого имеется аналитическое решение. Метод применяется к расчету вытекающих мод эллиптического ДВ, круглого и эллиптического ДВ типа «канал». Проводятся исследования практически важных волоконно-оптических устройств, работающих на вытекающих модах - трехслойных брэгговских волокон
Особенности программной реализации метода
Перечислим специфические особенности программной реализации метода для случая вытекающих мод по сравнению с вариантом метода, описанным в главе 1. Задание геометрии системы и триангуляция Для случая направляемых мод для задания геометрии системы была разработана программа - препроцессор, для решения задачи триангуляции была использована программа Tnangle (см.разделы 1.6.1, 1.6.2). В случае вытекающих мод для учета комплексного характера величины є(х,у) препроцессор был модифицирован, изменений программы Tnangle не требуется. Формирование матриц и решение проблемы собственных значений Для случая направляемых мод для решения этих задач была разработана программа - решатель СОММА (см.раздел 1.6.3). Рассмотрим модификации, внесенные в программу СОММА для расчета вытекающих мод. Вычисления производятся в комплексных величинах. Соответственно, используемые библиотеки ARPACK и UMFPACK также были заменены на версии с комплексной арифметикой. Проблема вычисления модифицированной функции Бесселя, обсуждаемая в разделе 2.2.2, была решена применением определенной техники интегрирования. Однако, при использовании описанной техники интегрирования требуется много машинного времени на расчет интегралов вида (2.33). Поэтому, описанная техника интегрирования была реализована на языке программы Mathematica 4.2, и использована в качестве эталона. В общедоступной библиотеке программного обеспечения NETLIB [79] был получен ряд подпрограмм, предназначенных для расчета модифицированной функции Бесселя второго рода для комплексных значений аргумента. Подпрограммы были протестированы, результаты расчета сравнивались с "эталонными" расчетами с использованием приведенной выше техники интегрирования. В результате этого отбора была выделена подпрограмма ZBESK, входящая в пакет AMOS [103 - 105]. Установлено, что комплексные значения, рассчитывамые с использованием этой подпрограммы, соответствуют как раз необходимому выбору ветви функции Макдональда.
В дальнейшем в программе-решателе использовалась подпрограмма ZBESK, что привело к кардинальному сокращению времени расчетов. Вывод и оценка результатов Для вывода и оценки результатов использовалась та же программа-постпроцессор, что и для случая направляемых мод (см.раздел 1.6.4). Была проведена модификация программы -постпроцессора для решения специфической проблемы, возникающей при разделении вещественной и мнимой компонент поля вытекающей моды. Действительно, собственные вектора обобщенной алгебраической задачи на собственные значения (1.20) определены с точностью до константы, в случае вытекающих мод — комплексной. Однако, в этом случае неизвестная комплексная константа, Е(х,у) - распределение напряженности поля вытекающей скалярной моды, С-Е(х,у) - величина, рассчитываемая программой - решателем. Это означает, что в полученном собственном векторе вещественная и мнимая компоненты поля представлены в виде линейной комбинации с неизвестными коэффициентами. Распределения вещественной и мнимой частей этой величины не имеют физического смысла. Поэтому для визуального представления поля вытекающей моды строилась величина положительной вещественной константы соответствует распределению средней за период объемной плотности электрической энергии вытекающей моды [32]. Соответствующий алгоритм был включен в программу - постпроцессор. Разработанная программа применялась для расчета вытекающих мод ДВ. Для систем, допускающих решение методом дисперсионного уравнения, была проведена серия тестовых расчетов. Результаты расчетов сравнивались с данными, полученными из решения дисперсионного уравнения. Рассчитанные величины были представлены в безразмерной форме. Волновые числа w и Р были отнормированы на характерный размер волновода b: wb, fib. Также использовался безразмерный параметр частоты: В случае, когда это определение неприменимо (ДВ типа «канал», ДВ типа «трубка», используется модифицированное определение V (см.соотв. разделы). Параметры V соответствующие частотам отсечки, снабжены индексом; Г... В качестве характерного размера b для круглого волновода используется радиус сердцевины, для эллиптического волновода - малая полуось сердцевины, для волновода с прямоугольной сердцевиной - половина меньшей стороны. Определим формат волновода f -alb как отношение характерных размеров сечения вдоль взаимно ортогональных осей, а - большая полуось волновода с эллиптической сердцевиной или половина большей стороны для волновода с прямоугольной сердцевиной, Ъ малая полуось волновода с эллиптической сердцевиной или половина меньшей стороны для волновода с прямоугольной сердцевиной. Расчетная область разбивалась на треугольные элементы, в типичном случае конечноэлементная сетка насчитывала несколько тысяч треугольников. Граница сопряжения Г, как правило, имела вид правильного 32-угольника и разбита на 100-200 отрезков. Для повышения точности расчетов использовалась равномерная сетка, граница Г приближена к сердцевине волновода. Генерировалась сетка с углами треугольных элементов не меньшими, чем 20 градусов. Такая сетка считается качественной, так как наличие в сетке треугольников с малыми углами снижает точность расчетов методом конечных элементов. Рассчитывались такие параметры волновода, как дисперсионная зависимость, потери на излучение и распределение поля - для заданных вытекающих мод. Дисперсионная зависимость представлена в виде Rej3b(V), а также Rewb(V) и lmwb(V). Потери на излучение вытекающей моды исследуются в виде величин Im/? и Im/?/Re/9 Величина Im/? описывает затухание вытекающей моды вдоль оси z. Величина Im/?/Re/? определяет уменьшение напряженности поля на одной длине волны. Достоинством величины Tm/?/Re/? является то, что она не зависит от используемой нормировки волнового числа Р. Распределение поля вытекающей моды представлено в видеі х, )!2, что соответствует распределению средней (за период) плотности энергии электрического поля моды. На монохромных рисунках темные области обозначают высокую, светлые - низкую плотность энергии моды. На цветных рисунках красный цвет обозначает
