ЧИСЛЕННЫЕ методы РАСЧЕТА ФОТОННЫХ структур

Андросик Андрей Борисович

к.т.н., доцент, Государственный открытый Университет, г. Москва

Воробьев Сергей Андреевич

к.т.н., доцент, Государственный открытый Университет, г. Москва

Мировицкая Светлана Дмитриевна

к.т.н., доцент, Государственный открытый Университет, г. Москва

E-mail: vorsa_57@mail.ru

Исходная точка развития вычислительной фотоники – анализ геометрии или топологии структуры интегрально-оптических элементов, т.е. определение характеристик пассивных и активных структур, включенных в схему. Анализ топологии структуры связан с нахождением функции оптического поля, удовлетворяющей уравнениям Максвелла или производным от них. Требуется удовлетворить граничные и начальные условия, а так же учесть свойства среды и условия возбуждения поля в волноводной структуре.

Численные методы (рис.1) преобразуют непрерывные интегральные или дифференциальные уравнения Максвелла в дискретную форму. Цель всех численных методов состоит в нахождении приблизительных решений уравнений Максвелла (или уравнений, полученных из них), которые удовлетворяют заданным граничным и начальным условиям.

Рисунок 1. Численные методы

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

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

Рисунок 2. Численные методы решения задач интегральной фотоники на основе интегральных уравнений Максвелла

Рисунок 3. Численные методы решения задач интегральной фотоники на основе дифференциальных уравнений Максвелла

Метод моментов или метод граничного элемента. Базовый метод моментов оперирует импульсными функциями расширения и тестовыми функциями Дирака. Этот подход был модифицирован с целью непосредственного получения интегральных уравнений  электрического (EFIE) и магнитного (MFIE) полей (рис.4).

Величины E и H являются функциями поля, а J – функция источника (текущая интенсивность). Инверсные операторы, таким образом, включают функции Грина, зависящие от граничных условий и распределения показателя преломления материала. В большинстве случаев эти интегральные уравнения формулируются в частотной области, хотя существуют применения также и для временной области. Вместо электрического E и магнитного H полей можно  найти скалярные или векторные потенциалы.

Рисунок 4. Метод моментов

Приближение дискретного диполя. Метод позволяет вычислять рассеяние и поглощение электромагнитных волн объектами произвольной геометрии и состава (рис.5). Исследуемый объект описывается областью поляризованных точек (диполей). Наилучшие результаты достигаются в случае, если размеры объекта сравнимы с длиной волны (область Ми). Материалы должны иметь |n-1|<1 t_o 3, n - комплексный показатель преломления. Идея метода заключается в разбиении малых объектов на малые объемы, каждый из которых несет момент диполя. Моменты диполя локального электрического поля рассматриваются отдельно для каждого диполя. Общее поперечное сечение вычисляется, используя приближение мультиполя для анизотропных и поглощающих рассеивателей. Поляризация малых объемов связана с диэлектрической функцией аппроксимацией Клаузиуса-Мозотти.

Рисунок 5. Приближение дискретного диполя

Метод конечных разностей. Метод конечной разности (рис.6) является наиболее часто используемым численным методом. Он осуществляет дискретизацию поперечного сечения анализируемой структуры; метод подходит для моделирования диэлектрических волноводов с произвольным поперечным сечением, которые могут быть изготовлены из изотропных однородных, неоднородных, анизотропных или материалов с потерями.

Рисунок 6. Метод конечных разностей

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

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

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

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

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

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

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

Метод конечных разностей во временной области - один из наиболее широко используемых числовых методов техники (рис.7). Это универсальный, методологически простой и наглядный метод. Волна, распространяющаяся через волноводную структуру, определяется путем прямого интегрирования во временной области (уравнение Максвелла) в дискретной форме. Главный недостаток метода конечных разностей во временной области – большое время вычислений.

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

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

Рисунок 7. Метод конечных разностей во временной области

Рисунок 8. Метод конечного элемента

Появление ложных решений, вызванных не выполнением условия дивергенции (Н=0), является серьезной  погрешностью метода. Избежать ложных решений можно путем формулирования вариационных выражений в членах различных компонент поля. Существуют различные формулировки, среди которых наиболее широко используется формулировка Н – поля для моделирования оптических волноводов в связи с наиболее легким выражением граничных условий. Однако такой подход не исключает возможности появления ложных решений. Формулировка в членах поперечных компонент электрического Е и магнитного Н полей не приводит к образованию ложных решений, однако она приводит к загущенности матрицы. Устранение ложных решений также можно обеспечить путем введения в вариационное выражение компенсирующего члена. Величина такого члена определяется эвристически или путем проверки условия нулевой дивергенции для каждой полученной моды и исключением мод, не удовлетворяющих этому условию. Описанные выше методики являются только частично успешными, поскольку ложные решения могут возникнуть во всем частотном спектре и отделить их от физически существующих мод практически невозможно. Наиболее удобно использовать методику граничного элемента, который гасит ложные решения непосредственно на нулевой частоте. В этом подходе, интерполяционные функции определяются как вектора, и удовлетворяется непрерывность тангенциальных компонент вдоль элементов. Непрерывность нормальных компонент поля не удовлетворяется, что приводит к ненулевой дивергенции. Однако, все ложные решения сводятся к нулевой частоте и, поэтому, легко идентифицируемы.

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

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

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

Для связи пространственной и спектральной областей было использовано быстрое преобразование Фурье (БПФ); такой метод получил название БПФ-МРЛ. Впоследствии метод БПФ-МРЛ был развит на случай слабо направляющих структур без учета векторных свойств поля. Использование параксиального приближения ограничивает метод структурами, в которых луч распространяется только в направлениях, составляющих малые углы  с осью распространения.

Чтобы избежать использования БПФ, параксиальное волновое уравнение было решено путем применения соответствующего вариационного выражения, использующего метод КЭ, а позже – с использованием метода КР, откуда и был получен метод конечной разности-распространения луча (КР-МРЛ). Сравнение методов БПФ – и КР-МРЛ показывает приемлемую точность при  большем шаге в случае КР. Также время расчета шага распространения в методе КР-МРЛ значительно меньше, что делает метод более привлекательным при исследовании комплексных структур.

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

Для моделирования условия открытой границы на конце области исследования, были применены поглощающие граничные условия (АВС). Впоследствии они были заменены  более эффективными условиями прозрачной границы (ТВС). Для поглощающих граничных условий нежелательные отражения абсорбируются материалом с потерями, расположенным на границе вычислительного окна. Наибольшим недостатком этой методики является то, что материал с поглощением вызывает определенные проблемы, поля на границе должны иметь нулевое значение; все это приводит к увеличению времени расчета и компьютерной памяти. В приближении прозрачной границы направление волны должно иметь определенный угол падения. Последний метод имеет меньше проблем в процессе реализации, более простой и не требует больших ресурсов памяти. Однако, размер окна должен все еще быть достаточно большим, чтобы не вызывать ослабления мощности части поля, которое распространяется в зоне сердцевины. Кроме того, эффективность методики прозрачной границы ограничена для структур с широким углом распространяющихся волн.

Рисунок 9. Метод распространения луча

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

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

1.            Андросик А.Б., Воробьев С.А., Мировицкая С.Д. Волноводная и интегральная  фотоника.- М.: МГОУ, 2011. 370 с.

2.            Андросик А.Б., Воробьев С.А., Мировицкая С.Д. Основы волноводной фотоники.- М.: МГОУ, 2009. 246 с.

3.            Андросик А.Б., Воробьев С.А., Мировицкая С.Д. Математические основы волноводной фотоники.- М.: МГОУ,  2010. 224 с.