Доклады Российской академии наук. Математика, информатика, процессы управления, 2020, T. 492, № 1, стр. 15-19

Бикомпактная разностная схема для уравнений Максвелла в слоистых средах

А. А. Белов 12*, Ж. О. Домбровская 1**

1 Московский государственный университет имени М.В. Ломоносова
Москва, Россия

2 Российский университет дружбы народов
Москва, Россия

* E-mail: aa.belov@physics.msu.ru
** E-mail: dombrovskaya@physics.msu.ru

Поступила в редакцию 27.09.2019
После доработки 27.09.2019
Принята к публикации 24.01.2020

Полный текст (PDF)

Аннотация

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

Ключевые слова: уравнения Максвелла, бикомпактные схемы, слоистые среды, условия сопряжения, дисперсия вещества

1. ПРОБЛЕМА

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

Многие известные методы (например матричные методы, методы конечных элементов, методы дискретных источников) дают хорошие результаты для достаточно узких классов задач. Наиболее универсальны разностные методы. Для них закрепилось общее название Finite-Difference Time Domain (FDTD) method [1]. Практически все они сводятся к схеме типа “крест” [2], в которой электрическое и магнитное поля относятся к целым и полуцелым пространственным узлам и временным слоям. Методам FDTD посвящено более 2000 работ, но принципиальные трудности этих методов не устранены.

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

Во-вторых, показатель преломления актуальных материалов (Si, Ge и др.) зависит от частоты сложным образом. Методы FDTD позволяют учитывать только простейшие модельные законы дисперсии (полиномиальный, модели Друде и Лоренца) [3]. Реальные законы дисперсии приближенно описывают комбинацией модельных законов. Погрешность такого приближения также трудно оценить. В результате нередко расчеты по стандартным пакетам методов FDTD плохо согласуются с экспериментами.

В данной работе на основе идеи бикомпактности [4, 5] построена консервативная разностная схема для системы одномерных уравнений Максвелла. Шаблон схемы включает только два узла по пространству, причем границы раздела выбираются узлами сетки. Такая схема обеспечивает второй порядок точности на решениях со слабыми и сильными разрывами. Предложен принципиально новый способ учета дисперсии среды. Закон дисперсии может быть произвольным. Построены физически содержательные и сложные для расчета примеры с точными решениями. Расчеты показали, что предложенная схема кардинально превосходит известные по точности и надежности.

2. РАЗНОСТНАЯ СХЕМА

Стационарная задача

1. Рассмотрим распространение плоской электромагнитной волны частоты ω перпендикулярно набору из Q плоско-параллельных пластин общей толщиной a. Пусть все объемные и поверхностные токи зависят от времени как $ \sim {\kern 1pt} {{e}^{{ - i\omega t}}}$.

Ориентируем координатную ось z перпендикулярно пластинам. Тогда в одномерной задаче электрическое поле ${\mathbf{E}}$ направлено по оси x, магнитное поле ${\mathbf{H}}$ – по оси y, волновой вектор ${\mathbf{k}}$ – по оси z. Здесь $k = \frac{\omega }{c}$, где c – скорость света. Диэлектрическая проницаемость ${{\varepsilon }_{q}}$, магнитная восприимчивость ${{\mu }_{q}}$, объемная плотность токов ${{{\mathbf{J}}}_{q}}$ $q$-й пластины, а также поверхностная плотность токов ${{{\mathbf{j}}}_{q}}$ на q-й границе не зависят от x и y. Вектора ${{{\mathbf{J}}}_{q}}$ и ${{{\mathbf{j}}}_{q}}$ направлены по оси x. Амплитуды полей зависят только от z. Одномерная постановка актуальна во многих приложениях (например, в интегральной фотонике).

2. Задача описывается интегральными уравнениями Максвелла

(1)
$\begin{gathered} \int\limits_\Gamma {{{{\mathbf{H}}}_{q}}} d{\mathbf{l}} = 4\pi {{c}^{{ - 1}}}\int\limits_S {{{{\mathbf{J}}}_{q}}} d{\mathbf{s}} + {{c}^{{ - 1}}}( - i\omega )\int\limits_S {{{{\mathbf{D}}}_{q}}} d{\mathbf{s}}, \\ {{{\mathbf{D}}}_{q}} = {{\varepsilon }_{q}}{{{\mathbf{E}}}_{q}}, \\ \end{gathered} $
(2)
$\int\limits_\Gamma {{{{\mathbf{E}}}_{q}}} d{\mathbf{l}} = - {{c}^{{ - 1}}}( - i\omega )\int\limits_S {{{{\mathbf{B}}}_{q}}} d{\mathbf{s}},\quad {{{\mathbf{B}}}_{q}} = {{\mu }_{k}}{{{\mathbf{H}}}_{q}}.$

Здесь $S$ – произвольная поверхность, ограниченная контуром $\Gamma $. Используется система единиц СГС. На внешних границах поставим условия излучения

(3)
$\begin{gathered} {{\partial }_{z}}{{{\mathbf{E}}}_{1}} + i\omega {{c}^{{ - 1}}}{{{\mathbf{E}}}_{1}} = 2ik{{{\mathbf{E}}}^{0}},\quad z = 0; \\ {{\partial }_{z}}{{{\mathbf{E}}}_{Q}} - i\omega {{c}^{{ - 1}}}{{{\mathbf{E}}}_{Q}} = 2ik{{{\mathbf{E}}}^{a}}{{e}^{{ - ika}}},\quad z = a. \\ \end{gathered} $

Здесь ${{{\mathbf{E}}}^{0}}$, ${{{\mathbf{E}}}^{a}}$ – заданные амплитуды волн, падающих справа и слева. Условия (3) аналогичны парциальным условиям излучения Свешникова [6], естественно учитывающим падение и отражение лазерного луча. По-видимому, такие условия используются впервые. Ранее падающую волну задавали моделью фиктивного источника. На границах слоев поставим условия сопряжения

(4)
$\begin{gathered} {{{\mathbf{e}}}_{z}} \times ({{{\mathbf{E}}}_{q}} - {{{\mathbf{E}}}_{{q - 1}}}) = 0,\quad {{{\mathbf{e}}}_{z}}({{{\mathbf{D}}}_{q}} - {{{\mathbf{D}}}_{{q - 1}}}) = 0, \\ {{{\mathbf{e}}}_{z}} \times ({{{\mathbf{H}}}_{q}} - {{{\mathbf{H}}}_{{q - 1}}}) = 4\pi {{c}^{{ - 1}}}{{{\mathbf{j}}}_{q}},\quad {{{\mathbf{e}}}_{z}}({{{\mathbf{B}}}_{q}} - {{{\mathbf{B}}}_{{q - 1}}}) = 0. \\ \end{gathered} $

3. Введем сетку $0 = {{z}_{0}} < {{z}_{1}} < \; \ldots \; < {{z}_{N}} = a$, ${{z}_{{n + 1}}}$ – – zn = Δzn так, чтобы границы слоев были узлами. Такие сетки называют специальными [7].

Для (1) в качестве поверхности $S$ выберем на плоскости y = 0 прямоугольник ${{z}_{n}}\,\, \leqslant \,\,z\,\, \leqslant \,\,{{z}_{{n + 1}}}$, $0\,\, \leqslant \,\,x\,\, \leqslant \,\,1$. Вычислим все интегралы в (1) по формуле трапеций. Это эквивалентно линейной интерполяции полей между узлами

(5)
$\begin{gathered} E(z) = {{\alpha }_{n}}z + {{\beta }_{n}},\quad H(z) = {{\gamma }_{n}}z + {{\delta }_{n}}, \\ {{z}_{n}} \leqslant z \leqslant {{z}_{{n + 1}}},\quad 0 \leqslant n \leqslant N - 1. \\ \end{gathered} $

Здесь ${{\alpha }_{n}}$, ${{\beta }_{n}}$, ${{\gamma }_{n}}$, ${{\delta }_{n}}$ – коэффициенты интерполяции. Сохраняя члены $O(\Delta {{z}^{2}})$ и сокращая $\Delta {{z}_{n}}$ в правой и левой частях равенства, получим

(6)
$\begin{gathered} {{\gamma }_{n}} = 4\pi {{c}^{{ - 1}}}{{J}_{q}}({{{\bar {z}}}_{n}}) - i\omega {{c}^{{ - 1}}}{{\varepsilon }_{q}}({{{\bar {z}}}_{n}})({{\alpha }_{n}}{{{\bar {z}}}_{n}} + {{\beta }_{n}}), \\ {{{\bar {z}}}_{n}} = 0.5\,({{z}_{n}} + {{z}_{{n + 1}}}). \\ \end{gathered} $

Для уравнения (2) в качестве поверхности $S$ выберем на плоскости x = 0 прямоугольник ${{z}_{n}} \leqslant z \leqslant {{z}_{{n + 1}}}$, $0 \leqslant y \leqslant 1$. Аналогично получим

(7)
${{\alpha }_{n}} = i\omega {{c}^{{ - 1}}}{{\mu }_{q}}({{\bar {z}}_{{}}}{{_{n}}_{n}})({{\gamma }_{n}}{{\bar {z}}_{n}} + {{\delta }_{n}}).$

4. Запишем условия сопряжения в каждом внутреннем узле, подставляя (5),

(8)
$\begin{array}{*{20}{c}} {{{\alpha }_{n}}{{z}_{{n + 1}}} + {{\beta }_{n}} = {{\alpha }_{{n + 1}}}{{z}_{{n + 1}}} + {{\beta }_{{n + 1}}},} \\ {{{\gamma }_{{n + 1}}}{{z}_{{n + 1}}} + {{\delta }_{{n + 1}}} - {{\gamma }_{n}}{{z}_{{n + 1}}} - {{\delta }_{n}} = 4\pi {{c}^{{ - 1}}}j({{z}_{{n + 1}}}).} \end{array}$

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

Остается аппроксимировать условия излучения (3). Выразим производную ${{\partial }_{z}}{{E}_{q}}$ из дифференциального уравнения Максвелла rotEq = $i\omega {{c}^{{ - 1}}}{{\mu }_{q}}{{{\mathbf{H}}}_{q}}$ и подставим интерполяцию (5) для H:

(9)
$\begin{array}{*{20}{c}} {i\mu ({{{\bar {z}}}_{0}}){{\delta }_{0}} + i\sqrt {\varepsilon ({{{\bar {z}}}_{0}})\mu ({{{\bar {z}}}_{0}})} {{\beta }_{0}} = 2i\sqrt {\varepsilon ({{{\bar {z}}}_{0}})\mu ({{{\bar {z}}}_{0}})} {{E}^{0}};} \\ \begin{gathered} i\mu ({{{\bar {z}}}_{N}})({{\gamma }_{N}}a + {{\delta }_{N}}) - \\ - \;i\sqrt {\varepsilon ({{{\bar {z}}}_{N}})\mu ({{{\bar {z}}}_{N}})} ({{\alpha }_{N}}a + {{\beta }_{N}}) = 2i\sqrt {\varepsilon ({{{\bar {z}}}_{N}})\mu ({{{\bar {z}}}_{N}})} {{E}^{a}}. \\ \end{gathered} \end{array}$

Равенства (6)–(9) представляют собой систему 4N алгебраических уравнений относительно 4N неизвестных ${\text{\{ }}{{\alpha }_{n}},{{\beta }_{n}},{{\gamma }_{n}},{{\delta }_{n}}{\text{\} }}$. Эта система является разностной схемой для стационарной задачи (1)–(4).

5. Схема выведена из законов сохранения, поэтому она консервативна. Принципиально новыми являются следующие аспекты.

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

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

Нестационарная задача

1. В постановке (1)–(4) следует заменить множитель $( - i\omega )$ на производную ${{\partial }_{t}}$. Объемные ${{J}_{q}}(z,t)$ и поверхностные ${{j}_{q}}(z,t)$ токи могут произвольно зависеть от времени. Слева и справа на рассеиватель падают волновые пакеты ${{E}^{0}}(t)exp( - i\omega t)$, ${{E}^{a}}\left( {t + \frac{a}{с}} \right){\text{exp}}\,( - i\omega t$ika) с несущей частотой $\omega $ и заданными огибающими ${{E}^{{0,a}}}$.

2. При распространении волнового пакета в линейной диспергирующей среде для каждой гармоники реализуются свои $\varepsilon $ и $\mu $ и своя скорость распространения. Поэтому пакет следует разбить на монохроматические компоненты, решить задачу для каждой из них в отдельности и просуммировать полученные решения. Разбиение пакета на компоненты сводится к прямому преобразованию Фурье, суммирование решений – к обратному.

3. Выполним численно преобразование Фурье падающих волновых пакетов, объемных и поверхностных токов по формуле трапеций, используя всюду одинаковые сетки по времени и наборы частот ${\text{\{ }}{{\omega }_{m}}{\text{\} }}$, ${{\omega }_{{m + 1}}} - {{\omega }_{m}} = \Delta {{\omega }_{m}}$. Для каждой частоты получим свои амплитуды компонент падающих волновых пакетов $E_{m}^{{0,a}}$ и источников поля ${{({{J}_{q}})}_{m}}$, ${{({{j}_{q}})}_{m}}$. Решим описанным выше методом набор стационарных задач, каждая из которых соответствует своей частоте ${{\omega }_{m}}$. Полученные амплитуды просуммируем по всем частотам ${{\omega }_{m}}$. Это и даст решение нестационарной задачи $E(z,t)$, $H(z,t)$.

Доказано, что описанный алгоритм сходится со вторым порядком точности. Для него также применимы ричардсоновские оценки точности и экстраполяционное уточнение. При этом сетки по $z$, $t$, $\omega $ должны сгущаться одновременно и в одинаковое число раз.

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

3. ТЕСТОВЫЕ РАСЧЕТЫ

1. Пусть граница раздела сред расположена в точке $z = b < a$. Слева от нее ${{\varepsilon }_{1}} = {{\mu }_{1}} = 1$, справа – ${{\varepsilon }_{2}} = 8$, ${{\mu }_{2}} = 2$. Падающие волны ${{E}^{{0,a}}}$ и объемные токи ${{J}_{{1,2}}}$ отсутствуют, а по границе раздела течет поверхностный ток ${{j}_{1}}$, излучающий волны в обоих направлениях оси z.

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

2. Точное решение этой задачи имеет вид

$\begin{gathered} E = - \xi {{j}_{1}}\left( {\frac{{[b - z]}}{{c - t}}} \right),\quad H = - E,\quad z < b, \\ E = \xi \sqrt {\frac{\varepsilon }{\mu }} {{j}_{1}}\left( {\frac{{[z - b]\sqrt {\varepsilon \mu } }}{c} - t} \right),\quad H = \sqrt {\frac{\varepsilon }{\mu }} E,\quad z > b. \\ \end{gathered} $

Здесь $\xi = 4\pi {{c}^{{ - 1}}}{{\left( {1 + \sqrt {\frac{\varepsilon }{\mu }} } \right)}^{{ - 1}}}$. Задача рассматривалась в двух постановках: стационарной при ${{j}_{1}} = {{e}^{{ - i\omega t}}}$ и нестационарной при ${{j}_{1}} = exp\left\{ { - \frac{{{{t}^{2}}}}{{t_{0}^{2}}} - i\omega t} \right\}$. В первом случае решение показано на рис. 1. Видно, что H имеет сильный разрыв, а E – слабый при z = b. Эта задача особенно трудна для численного расчета. В литературе она не рассматривалась.

Рис. 1.

Решение (10) в стационарном случае, вертикальная линия – граница раздела сред.

Погрешность определялась как разность численного решения с точными (10). Зависимость погрешности от числа шагов для разных схем приведена на рис. 2. Наклон графиков наших схем показывает сходимость со вторым порядком точности. Точность 1% достигается уже при N ~ 200 узлов сетки.

Рис. 2.

Погрешности; ⚪ – стационарная схема, △ – нестационарная схема, ◻ – явная схема “крест”.

3. Для сравнения был проведен расчет по явной схеме “крест” метода FDTD. Погрешность приведена на рис. 2. Она огромна ∼200% и не убывает при уменьшении шага; т.е. сходимости нет! Тем самым, общеупотребительная схема “крест” абсолютно непригодна для расчета слоистых сред.

4. Были проведены расчеты нестационарных задач в диспергирующих средах, показавшие аналогичные результаты.

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

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

  1. Inan U.S., Marshall R.A. Numerical electromagnetics. The FDTD method. Cambridge: Cambridge Univ. Press, 2011.

  2. Yee K.S. // IEEE Trans. Antennas. Propag. 1966. V. 14. № 3. P. 302.

  3. Sullivan D.M. Electromagnetic simulation using the FDTD method. IEEE Press, 2000.

  4. Калиткин Н.Н., Корякин П.В. // ДАН. 2007. Т. 419. № 6.

  5. Калиткин Н.Н., Корякин П.В. // Матем. моделирование. 2009. Т. 21. № 8. С. 44.

  6. Свешников А.Г. // ДАН. 1950. Т. 3. № 5. С. 517–520.

  7. Толстых А.И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. М.: Наука, 1990.

  8. Калиткин Н.Н., Альшин А.Б., Альшина Е.А., Рогов Б.В. Вычисления на квазиравномерных сетках. М.: Физматлит, 2005.

Дополнительные материалы отсутствуют.

Инструменты

Доклады Российской академии наук. Математика, информатика, процессы управления