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

08/01/2018

Полноволновое моделирование больших волновых оптических систем при строгом решении системы уравнений Максвелла классическими методами является очень ресурсоёмким. Это связано с требованием использования масштабной подробной и густой сетки, которая должна разрешать длину волны в среде. Технология Beam Envelopes (метод огибающих пучка), реализованная в программном обеспечении COMSOL Multiphysics®, позволяет сократить затраты вычислительных ресурсов для ряда задач волновой оптики. В данной статье нашего корпоративного блога, мы обсудим использование физического интерфейса Electromagnetic Waves, Beam Envelopes, его преимущества, ограничения и приложения к комплексным оптическим расчётам.

Сравнение методик расчета масштабных моделей волновой оптики

При электромагнитном моделировании для получения точного решения системы уравнений Максвелла длина волны всегда должна быть корректно разрешена конечно-элементной сеткой. Это требование затрудняет работу с моделями, которые являются большими по сравнению с длиной волны. Существует несколько альтернативных подходов к решению задач волновой оптики в частотной области, которые позволяют исследовать большие модели. Среди них можно выделить так называемые дифракционные формулы, в т.ч. формулы Фраунгофера, Френеля-Кирхгофа и Рэлея-Зоммерфельда, а также метод распространения пучка (beam propagation method) и его модификации, в т.ч. параксиальный метод распространения пучка и метод углового спектра (Ref. 1).

Большинство указанных выше альтернативных методов используют определенные приближения/упрощения для уравнения Гельмгольца. Они позволяют обрабатывать большие модели за счет так называемого метода распространения, для которого данные о поле в интересующей плоскости получаются на основе из известного поля в другой плоскости. Таким образом, вам не нужно делать дискретизацию для всей расчетной области (между указанными плоскостями), т.к. просто нужна 2D-сетка для интересующей плоскости.

По сравнению с этими техниками интерфейс Electromagnetic Waves, Beam Envelopes в COMSOL Multiphysics (далее будем коротко его называть интерфейсом Beam Envelopes) решает точное решение уравнения Гельмгольца в области без каких-либо аппроксимаций. Он позволяет обрабатывать большие модели, а требование к сетки может быть значительно ослаблено при выполнении определенных условий (или ограничений).
Линза, которая моделируется с помощью метода технологии Beam Envelopes.
Пучок с длиной волны 1 мкм, проходящий через линзу с фокусным расстоянием в миллиметровом диапазоне. Расчет с помощью технологии Beam Envelopes.

Ниже мы более подробно обсудим интерфейс Beam Envelopes.

Теория, лежащая в основе технологии Beam Envelopes

Давайте взглянем на математику, которая реализована в интерфейсе Beam Envelopes. Если вы добавите этот интерфейс в модель, выберете его корневой узел и измените Type of phase specification на User defined, то в разделе Equation вы увидите следующую формулу:

(\nabla-i \nabla \phi_1) \times \mu^{-1}_r (( \nabla-i \nabla \phi_1) \times
{\bf E1}) -k_0^2 \left (\epsilon_r -\frac{j \sigma}{\omega \epsilon_0} \right ) {\bf E1} = 0

Здесь \bf E1– это базовая переменная, для которой проводится расчет. Она как раз является огибающей.

В векторном представлении поля \bf E1 соответствует (медленно меняющейся) амплитуде, а \phi_1 – обобщенной фазе,

{\bf E} = {\bf E1} \exp(-i \phi_1).

Таким образом первое уравнение – управляющее уравнение для интерфейса Beam Envelopes – может быть получено путем подстановки второго определения для электрического поля в уравнение Гельмгольца. Если мы знаем \phi_1, то единственная неизвестная переменная – это \bf E1, и мы можем решить задачу. Подчеркнем, что закон изменения обобщенной фазы \phi_1 в пространстве должен быть задан a priori.

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

Распространенный вопрос: "Когда понадобятся именно данные об огибащей, а не о самом поле?" Одним из таких приложений является моделирование линз. В таких задачах часто требуются данные по интенсивности, а не комплексное электрическое поле. Фактически интенсивность – это квадрат нормы огибающей. Поэтому получение данных о распределении огибающей достаточно для последующего анализа.

Что произойдет, если обобщенная фазовая функция точно не известна?

Используемая формулировка технологии Beam Envelopes генерирует дополнительные вопросы

  • Что делать, если обобщенная фаза точно не известна?
  • Можно ли использовать интерфейс Beam Envelopes в таких случаях?
  • Насколько верными окажутся при этом результаты?

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

Одномерный случай

Рассмотрим простейший тестовый случай: плоская волна, Ez = \exp(-i k_0 x), где k_0 = 2\pi / \lambda_0 для длины волны \lambda_0 = 1 мкм, распространяется в прямоугольной области длиной 20 мкм. Мы намеренно используем короткий домен в иллюстративных целях.

Волна (с поляризацией вне плоскости) входит в область с левой границы и проходит через правую границу без отражения. Это можно описать в интерфейсе Beam Envelopes, добавив граничное условие Matched Boundary Condition с активированным возбуждением слева и без возбуждения справа, а также добавив граничное условие Perfect Magnetic Conductor сверху и снизу (то есть мы не заботимся о направлении y и считаем, что среда в этом направлении фактически бесконечна).

Правильные настройки при задании фазы показаны на рисунке ниже.
Скриншот графического интерфейса COMSOL Multiphysics, показывающий настройки раздела Wave Vectors интерфейса Beam Envelopes.

Мы, очевидно, получим ответ Ez = \exp(-i k_0 x), зная, что правильная фазовая функция равна k_0 x или волновой вектор a priori задан как (k_0,0). Подставляя фазовую функцию во второе уравнение, мы в свою очередь получаем E1z = 1, т.е. постоянную функцию.

Сколько элементов сетки нам потребуется для разрешения постоянной функции? Только один! См. заметку нашего блога о высокочастотном моделировании.

Приведенные ниже результаты показывают распределение в среде огибающей \bf E1 и нормы \bf E, ewbe.normE, которая равна |{\bf E1}|. Видно, что при задании корректного выражения для фазы мы получаем правильную огибающую функцию, постоянную единицу в данном случае, для любого числа сеточных элементов. В качестве дополнительной информации на графике отображается также фаза \bf E1z, arg(E1z). Она равна нулю, как и ожидалось.

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

Теперь давайте посмотрим, что произойдет, если наше предположение о фазовой функции будет немного ошибочным — допустим, (0.95k_0,0) вместо точного значения (k_0,0). Какое решение мы получим в этом случае? Давайте посмотрим:
Три результата (при разной сетке) для огибающей функции, нормы электрического поля и фазы огибающей функции при задании некорректного закона изменения фазы
Огибающая функция (красная линия), норма электрического поля (синяя линия) и фаза огибающей функции (зеленая линия) для некорректно заданной фазовой функции, 0.95 k0x при различной густоте сетки.

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

Мы знали, что верный ответ описывается как Ez = \exp(-i k_0 x), но мы "намеренно" использовали неверную оценку для фазовой функции в настройках интерфейса в программном обеспечении COMSOL®. Подставляя неправильную фазовую функцию во второе уравнение, мы получим \exp(-i k_0 x)= {\bf E1z} \exp(-0.95i k_0 x). Получается, что {bf E1z} = \exp(-0.05i k_0 x), т.е. она больше не является константой. Это волна с характерным масштабом \lambda_b= 2\pi/0.05k_0 = 20 мкм, которая часто называется длиной волны биения.

Давайте взглянем на график, приведенный выше, для сетки из шести элементов. Мы получаем именно то, что и ожидалось (красная линия), то есть {\bf E1z} = \exp(-0.05i k_0 x). На графике по умолчанию отображается действительная часть переменной, т.е. {bf E1z} = \cos(-0.05 k_0 x). Графики для менее подробного сеточного разрешения показывают лишь приближенное решение для огибающей функции. Это вполне ожидаемо для конечно-элементного моделирования: более грубая сетка дает менее точные и приближенные результаты.

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

А что насчет нормы поля \bf E? Посмотрите на синюю линию на графике выше. Похоже, что программное обеспечение COMSOL Multiphysics сгенерировало правильное решение для ewbe.normE, т.е. постоянное распределение для нашего случая. Давайте это проверим: подставляя во второе уравнение как неправильную (аналитическую) фазовую функцию, так и неправильную огибающую с биением, мы получаем {\bf Ez} = \exp(-0.05i k_0 x) \times \exp(-0.95i k_0 x) = \exp(-i k_0 x), что является правильным быстро-осциллирующим полем!

Если мы возьмем его норму \bf E, то получим правильное решение. Это иммено то, чего мы хотели добиться. Обратите внимание, что мы не можем корректно отобразить непосредственно \bf E, т.к. область может быть слишком большой (и не хватит сеточного разрешения), но мы можем найти \bf E аналитически и визуализировать норму \bf E даже с грубой сеткой.

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

Независимо от того, какую фазовую функцию мы используем в COMSOL Multiphysics, важно, чтобы правильно решалось первое уравнение для \bf E1а фазовая функция была непрерывна по области. Если в расчетной области находится несколько материалов, то непрерывность фазовой функции на границах их раздела также критична для получения точного решения. Возможно мы вернемся к этому в будущем, но исходные комментарии уже можно получить в одной из наших прошлых заметок о высокочастотном моделировании.

Двумерный случай

До сих пор мы обсуждали случай скалярного волнового числа. В более общем случае фазовая функция задается волновым вектором. Некорректно заданные приближения для волнового вектора (при использовании технологии Beam Envelopes) приведут к "векторным" последствиям. Предположим, что у нас в системе та же плоская волна из первого примера, но теперь мы сделаем неверное предположение о фазе, то есть зададим k_0(x \cos \theta + y \sin \theta) вместо k_0 x . В этом случае модуль волнового числа верный, но по его компонентам мы используем неверную оценку. На этот раз биение будет двумерным.

Давайте начнем с тех же аналитических оценок, что и для одномерного случая. Мы имеем \exp(-i k_0 x)= {\bf E1z}(x,y) \exp(-i k_0 (x \cos \theta+y \sin \theta) ), функция огибающей теперь вычисляется как {\bf E1z} (x,y) = \exp(-i k_0 (x (1-\cos \theta) -y \sin \theta) ) , т.е. которая представляет собой наклонную волну, распространяющуюся в направлении (1-\cos \theta, -\sin \theta), с волновым числом биения k_b = 2 k_0/\sin (\theta/2) и длиной волны биения \lambda_b=\lambda_0/(2\sin (\theta/2)).

На следующих графиках представлены результаты для θ = 15° и области с размерами 3,8637 мкм х 29,348 мкм при различной густоте сетки. Задаются те же граничные условия, что и в рассмотренном выше одномерном случае. Единственное отличие состоит в том, что падающая волна на левой границе задана как {\bf E1z} (0,y) = \exp(i k_0 y \sin \theta) . (Обратите внимание, что мы должны задать соответствующее неправильному исходному предположению о фазовой функции неправильное граничное условие.)

В результате для самой густой сетки (справа) мы можем видеть, что \bf E1z вычисляется точно так же, как мы предполагали на основе аналитической оценки, а норма \bf Ez получается на выходе постоянной. Эти результаты согласуются с одномерным примером.

Ряд результатов с данными о норме электрического поля и огибающей функции для неправильно заданной фазовой функции.
Норма электрического поля (вверху) и огибающая функция (внизу) для неправильной фазовой функции k_0(x \cos\theta +y \sin\theta ), вычисленные при различной густоте сетки. Радужная цветовая шкала соответствует значению поля от -1 до 1.

Моделирование линзы с использованием интерфейса Beam Envelopes

Наконец мы дошли до основной задачи нашей заметки, а именно до моделирования прохождения электромагнитного пучка через оптические линзы миллиметровых размеров с помощью интерфейса Beam Envelopes. Что для этого понадобится? Мы уже обсудили, как получить правильное решение. Постановка задачи для итогового примера следующая: плоский пучок с фиксированной аппертурой падает на на плосковыпуклую линзу с радиусом кривизны 500 мкм и показателем преломления 1.5 (фокусное расстояние линзы составляет приблизительно 1 мм).

В качестве приближения фазовой функции мы используем \phi_1 = k_0 x, что совсем не точно. В области перед линзой происходит отражение, которое приводит к интерференции. В линзе очевидно наблюдаются многочисленные переотражения. После прохождения линзы, фаза распределяется по сферическому закона, а пучок фокусируется в небольшое пятно. Таким образом, эта фазовая функция сильно отличается от того, что происходит в расчетной области. И все же у нас есть подсказка. Если мы построим график для \bf E1z, то увидим биение.

Расчет длины волны биения внутри линзы.
График распределения \bf E1z. Вставка с крупным планом показывает биения внутри линзы.

Отчетливо видно, что в линзе происходит заметное биение (см. вставку на графике выше).

Обратите внимание, что начиная с версии 5.5, в модель рассматриваемой линзы (которая доступна в Библиотеке моделей и приложений) были внесены изменения, в т.ч. было добавлено граничное условие Transition boundary condition, которое имитирует антиотражающее покрытие, которое устраняет биения в линзе из-за переотражений (а следовательно и необходимость сгущать сетку). См. более новую заметку нашего блога о полноволновом моделировании линз. При этом сгущение стеки все равно понадобится для случая ненормального падения, для которого уже не получится использовать антиотражающее покрытие. Кроме того, в окрестности границы линзы все равно будет наблюдаться небольшое биение.

На самом деле, самая малая по величине длина волны биения – это \lambda_0/2 для области перед линзой. Чтобы доказать это, мы можем выполнить те же вычисления, что и в ранее рассмотренных примерах. Самая малая длина волны биения обусловлена интерференцией между падающим пучком и отраженным фронтом, но мы можем это проигнорировать, т.к. указанный эффект не влияет на прямое распространение через линзу. Т.е. мы можем пренебречь тем фактом, что наша сетка не будет разрешать биения перед линзой.

Длина волны биения в линзе равна 3\lambda_0/2 для волны, идущей в обратном направлении (отрицательное направление горизонтальной оси), и 2\lambda_0 для прямой волны (положительное направление горизонтальной оси) при n = 1.5, что также несложно вывести аналитически. Опять же, мы проигнорируем идущий в обратном направлении сигнал. На графике ниже видны биения с периодом 2\lambda_0 для волны в прямом направлении. Волна в обратном направлении составляет лишь небольшую часть (примерно 4% от падающего пучка для n= 1.5), поэтому она не идентифицируется на графике. На следующем рисунке показана сетка, разрешающая биение внутри линзы с помощью 10 сеточных элементов.

Сетка для разрешения длины волны биения внутри линзы, созданная в COMSOL Multiphysics.
Биения внутри линзы, разрешаемые с помощью 10 сеточных элементов на длину волны биения.

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

Вычисленное распределение нормы поля в итоговой модели показано на самом верхнем изображении в данной статье. Чтобы верифицировать результат, мы можем сначала рассчитать поле на выходной поверхности линзы с помощью стандартного интерфейса Electromagnetic Waves, Frequency Domain, а затем на основе формулы дифракции Френеля вычислить поле в плоскости фокуса. Результаты по норме поля в фокусе очень хорошо согласуются.

1D-график со сравнением результатов полученных на основе технологии Beam Envelopes и формулы дифракции Френеля.
Сравнение результатов, полученных с помощью интерфейса Beam Envelopes и с помощью классического расчета, дополненного оценкой на основе формулы дифракции Френеля. В первом случае используется сетка, которая разрешает биения внутри линзы с помощью 10 элементов на длину волны биения.

На следующем графике приведена зависимость итогового результата от размера сетки внутри линзы. Мы получаем довольно хороший результат с нашей стандартной рекомендацией \lambda_b/6, которая в данном примере сводится к \lambda_0/3. В любом случае это очень "экономичное" требование для сетки в области линзы.

1D-график, показывающий зависимость нормы поля на определенном расстоянии от линзы от заданного размера сетки.
Зависимость распределения нормы поля в фокусе от размера сетки внутри линзы.

Начиная с версии 5.3a в программном обеспечении COMSOL® также доступна учебная модель расчёта линзы Френеля, построенная с помощью интерфейса Beam Envelopes. Линзы Френеля обычно чрезвычайно тонкие (порядка длины волн). Даже при наличии дифракции на поверхности линзы и вокруг нее мелкая сетка вокруг части линзы существенно не повлияет на общее количество элементов сетки и степеней свободы.

Заключение

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

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

Дальнейшие шаги

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

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

  1. J. Goodman, Fourier Optics, Roberts and Company Publishers, 2005.

Комментарии (0)

Оставить комментарий
Войти | Регистрация
Загрузка...
РУБРИКАТОР БЛОГА COMSOL
РУБРИКИ
ТЕГИ