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

Полный текст

1. Введение Распространение электромагнитных волн в полых волноводах было описано в символьном виде в ставшей уже классической работе [1]. Исследование закрытых волноводов со сложным заполнением ϵ и m требует применения численных методов [2]. В работах [3, 4] был предложен метод исследования задачи волноводного распространения излучения, основанный на декомпозиции полей по четырём потенциалам. Этот подход может быть применён и к исследованию открытых волноведущих систем путём введения виртуальных стенок [5, 6]. В настоящей работе предложенный подход применяется к 1) расчёту нормальных мод для модельных задач и 2) расчёту коэффициентов отражения и прохождения волноводных мод на стыке двух волноводов. Результаты численных экспериментов верифицировались в работе путём сравнения результатов расчёта модельных примеров в рамках двух независимых реализаций метода четырёх потенциалов, первая реализация разработана в системе компьютерной алгебры Sage, вторая - в Maple. Данные, полученные в рамках численных экспериментов по расчёту прохождения излучения через стык волноводов, сравнивались с результатами расчётов аналогичной конфигурации методом поперечных сечений, полученных в работе Шевченко [7]. Вычисления проводились в системах компьютерной алгебры Maple, Sage: матричные коэффициенты задачи на собственные значения вычислялись символьно, задача на собственные значения решалась численно с использованием встроенных в систему компьютерной алгебры численных методов. 2. Описание поля в волноводе при помощи четырёх потенциалов Рассмотрим волновод постоянного односвязного сечения S с идеально проводящими стенками; относительно заполнения ε,m которого не будем делать пока никаких предположений. Ось Oz направим по оси цилиндра, нормаль к S будем обозначать как n®, касательный вектор, перпендикулярный к e®z как t®. Положим С = (¶x,¶y)T,Сў = (-¶ y,¶x)T. Декомпозиция Гельмгольца [4, теоремы 1 и 2] позволяет искать решение системы уравнений Максвелла в виде E®^ = Сue + Сўv e,H®^ = Сvh + Сўu h. (1) Введённые здесь четыре скалярные функции будем называть потенциалами и всегда предполагать, что они удовлетворяют граничным условиям ue = uh = n ЧСve = n ЧСvh¶S = 0 (2) на границе волновода. При этом автоматически выполняются граничные условия идеальной проводимости, а сами уравнения Максвелла можно записать в следующем виде [3, 4] Сў¶ zvh - 1 ikmDve + ikεСўv e -С¶zuh + ikεСue = 0, Сў¶ zue + 1 ikεDuh - ikmСўu h -С¶zve - ikmСvh = 0. (3) Если ε и m являются постоянными, эта система распадается на пару уравнений Гельмгольца, данный подход был преложен в работах [3, 4]. Мы воспользуемся некоторыми теоретическими результатами из этих работ, чтобы в данной статье привести результаты компьютерных экспериментов, содержащих символьные и численные расчёты в системе компьютерной алгебры Maple. 3. Дискретизация уравнений Максвелла Рассмотрим вначале однородный вдоль оси z волновод квадратного сечения S = {x О [0;1],y О [0;1]}. Пусть на некотором отрезке a < z < b заполнение ε,m не меняется вдоль оси. Дискретизацию по пространству будем проводить, используя стандартный базис метода Галёркина [8, 9]. Раскладывать приближенное решение будем по конечному числу функций, удовлетворяющих граничным условиям. В качестве базиса в данном случае подходит базис f®j j=14N2-2 , составленный из собственных функций оператора Лапласа с условиями Дирихле и условиями Неймана, а именно: f®j j=1N2 -1 = cos pnxcos pmy,0,0,0T,m,n = 0,†,N - 1,m + n > 0, (4) f®j j=N2-12N2-2 = 0,cos pnxcos pmy,0,0T,m,n = 0,†,N - 1,m + n > 0, (5) f®j j=2N2-23N2-2 = 0,0,sin pnxsin pmy,0T,m,n = 1,†,N, (6) f®j j=3N2-24N2-2 = 0,0,0,sin pnxsin pmyT,m,n = 1,†,N. (7) Приближенное решение системы (??) будем строить по аналогии с неполным методом Галёркина [2, 9] в виде разложения vh,ve,ue,uhT = е j=14N2-2 wj zf®j x,y. (8) После подстановки приближенного решения в уравнения системы и применения проекционной схемы метода Галёркина получим систему обыкновенных дифференциальных уравнений относительно вектора искомых коэффициентных функций w® = w1,w2,†,w4N2-2T: Bdw® dz + ikAw® + 1 ikCw® = 0®, (9) где A,B и C - квадратные матрицы, элементы которых представляют собой двойные интегралы по сечению волновода S: aij = в€Sε Сўf j2 + Сfj3 ЧСўf i1 -Сfi4 dxdy+ + в€Sm Сўf j4 + Сfj1 ЧСўf i2 -Сfi3 dxdy, (10) bij = в€SСўf j1 ЧСўf i1 + Сўf j3 ЧСўf i3 + Сfj2 ЧСfi2 + Сfj4 ЧСfi4 dxdy, (11) cij = в€SDfj2Dfi1 m -Dfj2Dfi3 ε dxdy, (12) где fjn, n = 1,†,4 - n-я компонента j-й вектор-функции f®j. 4. Волновод с прямоугольной вставкой Широкий класс волноводов - это волноводы с кусочно-постоянным заполнением ε и m. В этом случае интегралы (??)-(??), определяющие матричные элементы, можно вычислить в символьном виде. Если ε полином, то вычислить интеграл (??) в символьном виде невозможно в силу теоремы Лиувилля [10]. Рассмотрим структуру матриц A, B и C, полученных символьно методом четырёх потенциалов для квадратного волновода с произвольной прямоугольной вставкой. pict (a) Структура матрицы A при N = 4 pict (b) Структура матрицы B при N = 4 pict (c) Структура матрицы C при N = 4 Рис. 1. Структура матриц A, B, C На приведённом рис. ?? чёрными квадратами обозначены символьные выражения, тождественно не равные нулю, белыми квадратами - символьные нули. Структура матрицы A разреженная, причём ненулевые символьные элементы матрицы расположены в блочно-шахматном порядке. Матрица C обладает также разреженной структурой с блочно-шахматным расположением ненулевых элементов. Матрица B обладает более простой - диагональной структурой в рамках выбранного базиса. Вычисление матричных элементов осуществлялось символьно в Maple с помощью встроенной функции int. Далее полученные выражения использовались для численных расчётов. Этот приём позволил избежать многократного повторения трудоёмкой процедуры численного интегрирования. 5. Нормальные моды Однородная система обыкновенных дифференциальных уравнений (??) допускает решения, которые зависят от z как eikbz, такие решения называют нормальными модами волновода. Подставляя вид решения w® z = y®eikbz в (??) и сокращая ненулевые множители, получим задачу на собственные значения и собственные векторы квадратной матрицы K Ky® = by®, где матрица K выражается через матрицы A,B и C следующим образом: K = B-1 A - 1 k2C. Величина b, входящая в степень экспоненты, называется коэффициентом фазового замедления моды, он определяет характер распространения моды. Моды с вещественным b > 0 представляют собой волны, бегущие вдоль оси волновода в положительном направлении, моды с b < 0 - в отрицательном направлении, их называют распространяющимися. Моды с мнимым b экспоненциально убывают в одном из направлений, их называют эванесцентыми. Вычисление собственных мод требует решения задачи на собственные значения и собственные векторы матрицы K. Благодаря диагональной структуре матрицы B в выбранном базисе матрицу K также можно вычислить символьно. Вид разреженной вычисленной символьно матрицы K аналогичен виду матриц A и C и приведён на рис. ??. Приведём далее результаты численных расчётов коэффициентов фазового замедления для квадратного волновода из материала с ε1 = 1 единичной ширины и высоты с квадратной вставкой px = py = p из материала с ε2 = 2 (см. рис. ??). pict Рис. 2. Структура матрицы K при N = 4 pict Рис. 3. Квадратный волновод с квадратной вставкой ширины p С увеличением толщины или изменением диэлектрической проницаемости вставки мнимые собственные значения трансформируются в действительные (см. [3]), что соответствует появлению новой распространяющейся моды. Проиллюстрируем описанное поведение собственных значений при изменении толщины диэлектрической вставки с p = 0,668 до p = 0,670, которые приведены на рис. ??. Ось абсцисс соответствует действительной части собственных значений, ось ординат - мнимой. Собственные значения bj, лежащие на оси абсцисс соответствуют направляемым модам, причём моды, бегущие в положительном направлении оси z, имеют bj > 0, а бегущие в отрицательном направлении - bj < 0. На мнимой оси расположены bj, отвечающие эванесцентным модам. При толщине p = 0,668 (левый график на рис. ??) в волноводе имеется три направляемые моды, бегущие в положительном направлении z, и три - бегущие в отрицательном. Квадратиком на рис. ?? обозначено мнимое собственное значение, которое при толщине p = 0,668 приближается к нулю и при толщине p = 0,670 (правый график на рис. ??) оно переходит на действительную ось, что соответствует появлению новой направляемой моды волновода. pict (a) p = 0,668 pict (b) p = 0,670 Рис. 4. Собственные значения Случаю p = 0 соответствует полый волновод, собственные значения для которого известны аналитически. В результате сравнения с численными результатами, полученными описанным методом, относительная погрешность при вычислении собственных значений в этом случае составила d < 2 Ч 10-10. Численный расчёт собственных значений производится средствами системы компьютерной алгебры Maple. С помощью функции ImportMatrix производится загрузка матрицы K, предварительно вычисленной символьно. Чтобы получить численные значения элементов матрицы, необходимо в этом случае в символьные выражения, зависящие от параметров вставки, подставить численные значения этих параметров, используя функцию subs. Для полученной матрицы K можно далее вычислить собственные значения, используя встроенную функцию Eigenvalues. Для матрицы K интересующей нас размерности 194 ґ 194 вычисления собственных значений и собственных векторов на современном компьютере с конфигурацией Intel(R) Core(TM) i7-6500U 2.50 ГГц, 8.00 ГБ ОЗУ занимают менее чем 0.5 сек. Описанный подход позволяет всесторонне исследовать спектральные характеристики волноведущих систем с кусочно-постоянным заполнением за приемлемое время. 6. Задача рассеяния на стыке двух волноводов Предложенный подход к представлению уравнений Максвелла через вспомогательные потенциалы с последующей дискретизацией позволяет достаточно просто в конструктивном виде сформулировать задачу волноводной дифракции. В качестве примера рассмотрим простейшую в рамках нашего подхода задачу из этого класса - задачу о дифракции волны на стыке двух волноводов. Пусть имеются два волновода, имеющие одинаковые сечения, но различное заполнение, и пусть они стыкуются при z = 0. Иными словами, пусть ε = ε1(x,y),z < 0, ε2(x,y),z > 0, m = m1(x,y),z < 0, m2(x,y),z > 0. Будем рассматривать случай дифракции волноводной моды, налетающей из левого (z < 0) волновода на правый (z > 0). Поле в левом и правом волноводах может быть представлено в виде суперпозиции нормальных мод, приближенное вычисление которых осуществляется методом четырёх потенциалов. При z < 0 задача отыскания нормальных мод редуцируется к задаче на собственные значение и собственные векторы матрицы K(1) K(1)y®(1) = b(1)y®(1), (13) а при z > 0 - к задаче K(2)y®(2) = b(2)y®(2). (14) Поле в левом волноводе состоит из падающей на стык моды с Re(b1(1)) > 0 и амплитудным коэффициентом A и отражённых мод с Re(bj(1)) < 0 и амплитудными коэффициентами отражения Rj, j = 1,†,2N2 - 1. Поле в правом волноводе состоит из прошедших мод с Re(bj(2)) > 0, характеризующихся амплитудными коэффициентами прохождения Tj, j = 1,†,2N2 - 1. Задача дифракции такой волны на стыке состоит в отыскании поля, равного w = Aw(1)(x,y)eikb1(1)z + е Re close="}" bj(1)<0Rjwj(1)(x,y)eikbj(1)z при z < 0 и w = е Re bj(2)>0Tjwj(2)(x,y)eikbj(2)z при z > 0, где wj(1)(x,y) и wj(2)(x,y) определяет поперечное распределение j-й нормальной моды левого и правого волноводов соответственно. Они определены как wj(1) x,y = е m=14N2-2 ymj(1)f® m x,y,wj(2) x,y = е m=14N2-2 ymj(2)f® m x,y, (15) где ymj(1) и ymj(2) - это m-я компонента j-го собственного вектора задачи (??) и (??) соответственно. На границе должно выполняться условие непрерывности тангенциальных составляющих полей: [E®^] = 0,[H®^] = 0. В силу единственности (с точностью до аддитивной константы) представления полей через потенциалы [3, 4], это условие выполняется, если потенциалы непрерывны. Таким образом, условие сопряжения полей ведёт к легко выписываемому условию на коэффициенты Rj,Tj Aw(1)(x,y) + е Re close="}" bj(1)<0Rjwj(1)(x,y) = е Re close="}" bj(2)>0Tjwj(2)(x,y). (16) Поскольку и wj(1)(x,y) и wj(2)(x,y) согласно (??) раскладываются по системе функций f®j j=14N2-2 для выполнения (??) коэффициенты при каждой функции f®j, j = 1,†,4N2 - 2, должны быть равны между собой. Приравнивая коэффициенты при каждой f®j, j = 1,†,4N2 - 2, в которые входят линейные комбинации Rj, Tj, получим систему из 4N2 - 2 линейных алгебраических уравнений относительно искомых коэффициентов отражения Rj, j = 1,†,2N2 - 1 и прохождения Tj, j = 1,†,2N2 - 1. Число неизвестных Rj,Tj в точности равно числу уравнений, поэтому полученная система имеет и притом единственное решение. Описанная процедура получения системы линейных алгебраических уравнений реализовывалась в системах компьютерной алгебры Maple и Sage для сравнения результатов расчётов модельного стыка, речь о котором пойдёт ниже. 7. Численное решение задачи дифракции на стыке закрытых волноводов Рассмотрим стык полого волновода ε1 = 1, m1 = 1 с волноводом градиентным, которому соответствует ε2 = 1 + dxy(1 - x)(1 - y), m2 = 1. Значение безразмерного параметра d определяет, насколько сильно различаются материалы левого и правого волноводов. Рассмотрим случай, когда из полого волновода на волновод градиентный налетает одна направляемая мода с b в‰ѓ 0,778 и единичной амплитудой A = 1. Будем численно оценивать её коэффициент отражения от стыка |R1(d)|2 для d О [1;26] при N = 2. PIC Рис. 5. Коэффициент отражения первой моды от стыка при N = 2, (а) результаты, полученные в работе В.В. Шевченко, (б) результаты, полученные в данной работе. Отражение при малых значениях параметра d незначительно. Увеличение параметра d влечёт большее отражение волноводной моды от стыка, что физически соответствует поведению моды на стыке волноводов из различных материалов. Сравнение реализаций символьно-численного метода расчёта коэффициентов отражения и прохождения в Maple и Sage для рассматриваемого случая показывает, что |R1Maple(1) - R1Sage(1)| < 0,002. Точность совпадения обусловлена различными реализациями символьно-численного метода, а также различной точностью численного расчёта интегралов в системах Maple и Sage. 8. Численное решение задачи дифракции на стыке открытых планарных волноводов Рассмотрим теперь второй пример, сделав предварительное замечание. Одним из первых методов исследования закрытых нерегулярных волноводов был разработан Б. З. Каценеленбаумом [11] - метод поперечных сечений, который был обобщён на открытые волноводы В. В. Шевченко [12]. Попытки обосновать метод поперечных сечений для открытых волноводов привели к серьёзным теоретическим трудностям, связанным с наличием непрерывного спектра задачи, описывающего излучательные моды открытого волновода. Волноводная задача для закрытого нерегулярного волновода тем временем имеет корректную постановку. В работе [5] было показано, что расчёт открытого волновода можно проводить в рамках модели так называемого «ящика Дирихле». Предложенная модель содержит предположение о возможности поместить открытый волновод в объемлющий полый закрытый волновод, границы которого удалены от реальных границ волноводного слоя и поставить корректную математическую задачу. В этом случае поведение направляемых и излучательных мод открытого волновода будет адекватно описываться в рамках предложенной приближенной модели, что подтверждается численными экспериментами, приведёнными в [5] и настоящей работе. В рамках модели «ящика Дирихле» будем проводить расчёт стыка открытых планарных волноводов, который подробно исследовался в работе В. В. Шевченко [7]. Рассмотрим аналогичную структуру, помещённую в ящик, и применим для расчёта предложенный в настоящей работе метод. В работе [7] рассматривается стык волноводов одинаковой толщины с различными диэлектрическими проницаемостями: из левого планарного волновода с ε0 = 2, m0 = 1 толщины d падает первая направляемая мода на правый планарный волновод ε1 = s Ч ε0, m1 = 1 такой же толщины, причём длина волны связана с толщиной соотношением d∕l = 0,4. Покровная среда для обоих волноводов - воздух ε = 1, m = 1. В рамках модели «ящика Дирихле» единичной ширины и высоты толщина волноводного слоя была принята d = 0,01, чтобы высота «ящика» была много больше толщины волноводного слоя. Ниже приведён график зависимости коэффициента прохождения направляемой моды в зависимости от безразмерного параметра s = ε1∕ε0 при N = 7 (рис. ?? (b)). pict (a) pict (b) Рис. 6. Коэффициент прохождения первой моды через стык при N = 7 Приведённая на графике зависимость демонстрирует качественную согласованность с результатами работы В. В. Шевченко [7] (рис. ?? (a)). Численные значения также имеют одинаковый порядок, что говорит об адекватности предложенного метода, в том числе для исследования открытых волноводов, помещённых в «ящик Дирихле». 9. Заключение Рассмотренные численные примеры показывают, что предложенная в первой части статьи схема сведения системы уравнений Максвелла в волноводе к системе двух связанных уравнений Гельмгольца, записанных в «гамильтоновой форме», приводит в рамках неполного метода Галёркина к очень удобной для численно-аналитического анализа системе обыкновенных дифференциальных уравнений. За рамками статьи осталось обсуждение возможностей, которые открывает указанный в статье произвол в выборе и упорядочении базиса. Можно полагать, что этот подход позволит создать простые в использовании для прикладных исследований комплексы программ для исследования волноводного распространения электромагнитных волн.

Об авторах

Анастасия Александровна Тютюнник

Автор, ответственный за переписку.
Email: tyutyunnik_aa@rudn.university

ассистент кафедры прикладной информатики и теории вероятностей РУДН

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