Радиотехника и электроника, 2023, T. 68, № 6, стр. 521-526

ВВЕДЕНИЕ

Задачи рассеяния электромагнитных волн на трехмерных анизотропных диэлектрических структурах имеют огромное значение как с теоретической, так и с практической точки зрения. Важной характеристикой задач подобного рода является соотношение длины волны электромагнитного поля и размеров рассеивающей структуры. Если отношение размеров структуры и длины волны достаточно велико, то используются асимптотические методы решения. Огромный вклад в развитие этих методов для широкого класса задач внесли Д.С. Лукин и его ученики – А.С. Крюковский, Е.А. Палкин, Д.В. Растягаев и др. [1–4]. Если длина волны соизмерима или больше размеров рассеивающего объекта, то необходима математически строгая постановка задачи. Исследование рассматриваемых задач имеет принципиальное значение в электродинамике.

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

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

К интегральным уравнениям Фредгольма 2‑го рода сводятся следующие классы задач: рассеяние на однородном изотропном диэлектрическом теле с гладкой границей (можно рассматривать также несколько гладких непересекающихся поверхностей раздела сред внутри тела) [6–8]; рассеяние в ограниченной изотропной среде, диэлектрическая проницаемость которой является всюду дифференцируемой функцией координат и нигде не обращается в ноль [7]. Используя теорию интегральных уравнений Фредгольма, были доказаны теоремы существования и единственности соответствующих задач. Однако в общей постановке, тем более при наличии анизотропии среды, трехмерные задачи рассеяния электромагнитных волн не могут быть сведены к классическим интегральным уравнениям Фредгольма.

Объемные сингулярные интегральные уравнения относительно неизвестного электромагнитного поля по области неоднородности Q появились в электродинамике в 70-е годы ХХ в. [9]. Эти уравнения описывают задачи электромагнитного рассеяния на трехмерных неоднородных и анизотропных диэлектрических структурах в самой общей постановке. Кроме того, уравнения могут использоваться при исследовании задач рассеяния, для которых классическая постановка либо невозможна, либо затруднена, например во фрактальных средах или диэлектрических телах с негладкой границей. Однако эти уравнения относятся к многомерным сингулярным интегральным уравнениям, теория которых была построена значительно позже и в меньших деталях по сравнению с теорией интегральных уравнений Фредгольма.

С использованием результатов С.Г. Михлина и его учеников по теории сингулярных уравнений [10, 11] нами впервые были исследованы объемные сингулярные интегральные уравнения электродинамики [5, 15, 16]: получены условия фредгольмовости уравнений, доказаны теоремы существования и единственности решения, исследованы классическая и обобщенная постановки задач. Математически строго получены условия существования резонансных анизотропных диэлектрических структур [14].

Отметим, что в работах Ю.Г. Смирнова и его учеников [12, 13] были рассмотрены на основе интегродифференциальных уравнений задачи рассеяния на трехмерных изотропных диэлектрических структурах и идеально проводящих поверхностях. С использованием теории псевдодифференциальных уравнений были изучены корректные постановки задач, доказаны теоремы существования и единственности.

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

где ${{\varepsilon }_{{nm}}}$ – компоненты тензора диэлектрической проницаемости $\hat {\varepsilon }$ в декартовой системе координат. Для изотропной среды условие (1) принимает вид $\varepsilon \left( x \right) \ne 0,x \in Q$.

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

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

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

В объемных сингулярных интегральных уравнениях параметры среды входят в ядра уравнений в простой форме, что делает возможным их численное решение. По-видимому, первое численное решение этих уравнений было приведено в работе [9]. Рассматривалась задача рассеяния плоской волны на неоднородном диэлектрическом параллелепипеде. Методом коллокации задача была аппроксимирована СЛАУ небольшой размерности, которая была решена прямым методом Гаусса. В дальнейшем уравнения использовались для численного решения рядом исследователей, однако широкого применения они пока не получили.

Опишем основные проблемы, которые возникают при численном решении этих уравнений. Для большинства реальных задач в силу трехмерности уравнений после дискретизации уравнений возникают СЛАУ очень большой размерности N, N $ \gg $ 1000. Очевидно, что использование прямых методов практически невозможно, поскольку в памяти компьютера необходимо хранить $M\sim {{N}^{2}}$ чисел, а для решения СЛАУ нужно выполнить $T\sim {{N}^{3}}$ арифметических операций. Поэтому возможно использование только итерационных методов. В этом случае параметры M и T оцениваются формулами

где ${{M}_{A}}$ – количество элементов массива чисел, которые необходимы для алгоритма умножения матрицы СЛАУ на вектор; L – количество итераций для получения решения с заданной точностью; ${{T}_{A}}$ – число арифметических операций для умножения матрицы СЛАУ на вектор. Для плотных матриц произвольного вида ${{M}_{A}} = {{N}^{2}}$, ${{T}_{A}}\sim {{N}^{2}}$. Однако для плотных матриц, обладающих определенными свойствами, эти величины могут быть значительно меньше.

В данной статье рассматривается метод численного решения объемных сингулярных интегральных уравнений. С использованием свойств ядер интегральных уравнений и алгоритмов быстрого дискретного преобразования Фурье конструируются алгоритмы со значениями ${{M}_{A}}$ и ${{T}_{A}}$, практически пропорциональными N [5, 17–19].

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим следующий класс задач. В трехмерной ограниченной области Q среда характеризуется кусочно-дифференцируемой тензор функцией диэлектрической проницаемости $\hat {\varepsilon }(x)$, а вне области диэлектрическая проницаемость изотропна и равна ${{\varepsilon }_{0}}$, и везде $\mu = {{\mu }_{0}}$, $\operatorname{Im} \,{{\varepsilon }_{0}} = 0,$ $\operatorname{Im} \,{{\mu }_{0}} = 0$, $\operatorname{Re} \,{{\varepsilon }_{0}} > 0,$ $\operatorname{Re} \,{{\mu }_{0}} > 0$. Требуется определить электромагнитное поле, возбуждаемое в данной среде внешним полем и создаваемое сторонним током ${{\vec {J}}^{0}}$ с временной зависимостью в виде множителя $\exp ( - i\omega \,t)$. Нашей целью является изучение решений уравнений Максвелла

удовлетворяющих условию излучения на бесконечности где ${{k}_{0}} = \omega \sqrt {{{\varepsilon }_{0}}{{\mu }_{0}}} $, $u$ – любая из декартовых компонент полей $\vec {E}$ или $\vec {H}$. Кроме того, на поверхностях разрыва параметров среды тангенциальные компоненты полей должны быть непрерывны.

Поставленная задача может быть сведена к объемному сингулярному интегральному уравнению относительно электрического поля в Q [5, 15]:

где G – функция Грина для уравнения Гельмгольца $R = \left| {x - y} \right|$, $x = ({{x}_{1}},{{x}_{2}},{{x}_{3}})$, $y = ({{y}_{1}},{{y}_{2}},{{y}_{3}})$, а p.v. обозначает сингулярный интеграл, т.е. интеграл за вычетом бесконечно малого шара в окрестности точек x и y.

2. РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЙ

В работах [5, 16] доказана следующая теорема единственности решения рассматриваемых задач.

Теорема 1. Пусть тензор функция $\hat {\varepsilon }(x)$ является кусочно-дифференцируемой функцией координат, а поверхности разрыва параметров среды – кусочно-гладкие. Далее, тензор-функция $\hat {\varepsilon }(x)$ удовлетворяет условию (1), а эрмитова тензор-функция ${{\left( {{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varepsilon } }}_{{}}}(x) - \hat {\varepsilon }{\text{*}}(x)} \right)} \mathord{\left/ {\vphantom {{\left( {{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varepsilon } }}_{{}}}(x) - \hat {\varepsilon }{\text{*}}(x)} \right)} {(2i)}}} \right. \kern-0em} {(2i)}}$ неотрицательно определена в каждой точке области Q. Тогда решение однородных уравнений Максвелла (3), удовлетворяющее условиям непрерывности тангенциальных компонент полей на поверхностях разрыва параметров среды и условию излучения на бесконечности, может быть только тривиальным.

Приведенное в теореме условие неотрицательности тензора ${{\left( {{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varepsilon } }}_{{}}}(x) - \hat {\varepsilon }{\text{*}}(x)} \right)} \mathord{\left/ {\vphantom {{\left( {{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varepsilon } }}_{{}}}(x) - \hat {\varepsilon }{\text{*}}(x)} \right)} {(2i)}}} \right. \kern-0em} {(2i)}}$ физически означает, что в среде нет генерации энергии, т.е. среда является пассивной. В изотропном случае условие принимает вид $\operatorname{Im} \,\varepsilon (x) \geqslant 0.$

Для ответа на вопрос о существовании решения рассматриваемых задач сначала необходимо выбрать подходящее для анализа функциональное пространство. Интегралы от квадрата модуля характеристик электромагнитного поля присутствуют в законе сохранения энергии. Значит, можно полагать, что пространство интегрируемых с квадратом вектор-функций является наиболее “физичным” для исследования интегральных уравнений задач электромагнитного рассеяния. Далее будем использовать гильбертово пространство ${{\vec {L}}_{2}}$ со скалярным произведением, определяемым формулой

Дадим несколько определений, используемых далее.

Определение. Пусть A – линейный оператор, действующий в гильбертовом пространстве H. Тогда оператор $A{\text{*}}$, который также определен в H, называется сопряженным к A, если следующее равенство

выполняется для всех $f,g \in H$.

Решения однородного уравнения $Au = 0$ будем называть нулями оператора A. Обозначим размерность подпространства нулей через $n{\kern 1pt} {\kern 1pt} {\kern 1pt} (A)$. Тогда $n\,(A*)$ – размерность подпространства нулей сопряженного оператора $A{\text{*}}$. Разность

называется индексом оператора A.

Определение. Линейный оператор A, действующий в гильбертовом пространстве, называется нетеровым оператором, если значения $n{\kern 1pt} {\kern 1pt} {\kern 1pt} (A)$ и $n\,(A*)$ конечны.

Определение. Линейный оператор A называется фредгольмовым оператором, если он нетеров и его индекс равен нулю.

Доказано следующее утверждение [5, 15].

Теорема 2. Пусть компоненты тензора $\hat {\varepsilon }(x)$ всюду гельдер-непрерывные функции координат. Тогда, для того чтобы оператор сингулярного интегрального уравнения (5) был нетеровым оператором в ${{\vec {L}}_{2}}(Q)$, необходимо и достаточно выполнение условия (1).

Для изотропной среды условие (1) принимает вид

Теорема 2 имеет принципиальное значение, поскольку с ее помощью доказываются многие важные утверждения.

Доказана следующая теорема существования и единственности решения [5, 16].

Теорема 3. Пусть выполняются условия теоремы 1. Кроме того, тензор-функция ${\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {{{{\hat {\varepsilon }}}_{r}}{\kern 1pt} (x) - \hat {I}} \right)$ может не иметь обратную функцию в Q только на множестве точек меры ноль, например на границе Q. Тогда существует и единственно решение сингулярного интегрального уравнения (5) в гильбертовом пространстве ${{\vec {L}}_{2}}(Q)$.

Из теоремы 3, следует, что практически для всех ограниченных пассивных диэлектрических структур Q, удовлетворяющих условию (1), решение задачи рассеяния существует и единственно.

Представляют интерес задачи рассеяния, при которых классическая постановка либо затруднена, либо невозможна. Доказано следующее утверждение при минимальных ограничениях на параметры среды [15].

Теорема 4. Пусть тензор-функция $\hat {\varepsilon }(x)$ – ограниченная функция координат в области Q, а эрмитов тензор ${{\left( {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varepsilon } (x) - \hat {\varepsilon }*(x)} \right)} \mathord{\left/ {\vphantom {{\left( {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varepsilon } (x) - \hat {\varepsilon }*(x)} \right)} {(2i)}}} \right. \kern-0em} {(2i)}}$ положительно определен в каждой точке Q. Тогда существует и единственно решение сингулярного интегрального уравнения (5) в гильбертовом пространстве ${{\vec {L}}_{2}}(Q)$.

Большой интерес представляют среды, в которых условие (1) не выполняется. В работе [14] доказано, что при невыполнении условия (1) в среде без потерь в области неоднородности Q, могут существовать ненулевые решения однородной задачи. Если среда анизотропная, то в области Q может находиться электромагнитная энергия и при этом нет излучения в окружающее пространство.

3. МЕТОД РЕШЕНИЯ РАССМАТРИВАЕМЫХ УРАВНЕНИЙ

Рассмотрим применение метода коллокации для численного решения объемного сингулярного интегрального уравнения (5). Отметим важное обстоятельство: ядра интегральных операторов уравнения (5) зависят только от разности декартовых координат точек x и y. Поэтому при дискретизации целесообразно это учитывать с целью получения матрицы СЛАУ, которая будет обладать соответствующими свойствами симметрии.

Для этого в прямоугольной декартовой системе координат введем конечное множество точек – сетку, так чтобы область Q целиком находилась в прямоугольном параллелепипеде $\Pi $ со сторонами ${{N}_{l}}\,{{h}_{l}}$, l = 1, 2, 3, где ${{h}_{l}}$ – шаги сетки по декартовым координатам. При этом параллелепипед $\Pi $ разбивается данной сеткой на ячейки (элементарные параллелепипеды) $\Pi (p)$, $p = ({{p}_{1}},\,{{p}_{2}},{{p}_{3}})$, ${{p}_{l}} = 0,...,{{N}_{l}} - 1$. Мы определяем область $\tilde {Q}$ как объединение ${{N}_{Q}}$ ячеек, центры которых лежат внутри области Q. Очевидно, что ${{N}_{Q}} \leqslant {{N}_{1}}{{N}_{2}}{{N}_{3}}$. Тогда выбор сетки определяется формой области Q и требованиями: искомое решение, внешнее поле и параметры среды слабо изменялись внутри ячеек. Узловые точки, в которых определяются значения искомых функций, будем задавать в центрах ячеек x(p) с координатами

Интегральное уравнение (5) может быть преобразовано к следующему виду [18]: Здесь $\Pi (\varepsilon )$ – параллелепипед, диаметр которого $\varepsilon $, с центром в точке x = y, а $\hat {\gamma }$ – диагональный тензор, зависящий от формы параллелепипеда [18]. Если параллелепипед – куб, то коэффициенты тензора $\hat {\gamma }$ равны, как и для шара, 1/3.

Интегральное уравнение (9) будем аппроксимировать СЛАУ с размерностью ~${{N}_{Q}}$ относительно значений неизвестного поля u(p) в узловых точках области Q

где $u(p) \equiv u(x(p));$ ${{u}_{0}}(p) \equiv {{u}_{0}}(x(p));$ $u(q) \equiv u(x(q))$.

Поскольку узловые точки находятся в центре ячеек, точность аппроксимации интегрального оператора ~h², где h вычисляется по формуле

В (10) $u(p)$ и ${{u}_{0}}(p)$ – векторы, а B(p–q) и $\eta (q)$ – матрицы с размерностью три, которые определяются согласно (9) формулами Здесь а ${{\alpha }_{n}}$ задаются формулами Матрица системы линейных уравнений (10), получившаяся после дискретизации интегрального уравнения, имеет форму блочно-теплицевой матрицы. Очевидно, что основные вычислительные затраты при умножении матрицы СЛАУ на вектор, связаны с вычислением сумм вида где трехмерные матрицы B(p–q) определяются формулами (12)–(14).

Для вычисления W(p) в узловых точках $x(p) \in Q$ требуется выполнить ~$N_{Q}^{2}$ арифметических операций. Для уменьшения числа арифметических операций будем применять технику быстрого умножения теплицевых матриц на вектор, основанную на быстром дискретном преобразовании Фурье [20].

Доопределим величины $V(q)$ нулем в точках x(q) параллелепипеда П, не принадлежащих области Q. Тогда (15) можно записать в виде

Очевидно, что при $x(p) \in Q$ значения W(p) из (15) и (16) совпадают. В (16) матричная функция дискретного аргумента $B(p)$ определена для значений Обозначим через ${{\Pi }_{2}}$ параллелепипед со сторонами $2{{N}_{1}}{{h}_{1}}$, $2{{N}_{2}}{{h}_{2}}$ и $2{{N}_{3}}{{h}_{3}}$. Продолжим матричную функцию дискретного аргумента $B({{p}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{p}_{3}})$ на все целочисленные значения ${{p}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{p}_{3}}$, полагая ее периодической по каждой переменной с периодами соответственно $2{{N}_{1}},{\kern 1pt} 2{{N}_{2}},2{{N}_{3}}$. При этом доопределим функции $B({{p}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{p}_{3}})$ нулем в точках $({{N}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{p}_{3}})$,$({{p}_{1}},{\kern 1pt} {{N}_{2}},{\kern 1pt} {{p}_{3}})$, $({{p}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{N}_{3}})$, где ${{p}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{p}_{3}}$ – любые целые числа. Далее доопределим вектор-функцию дискретного аргумента $V({{p}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{p}_{3}})$ нулем во всех узловых точках ${{\Pi }_{2}}$, не принадлежащих $\Pi $, и продолжим ее на все целочисленные значения ${{p}_{1}},{\kern 1pt} {{p}_{2}},{\kern 1pt} {{p}_{3}}$, полагая ее периодической по каждой переменной с периодами соответственно $2{{N}_{1}},{\kern 1pt} 2{{N}_{2}},2{{N}_{3}}$.

Запишем выражение

С учетом изложенного ясно, что при $x(p) \in \Pi $ функция $W({{p}_{1}},\,{{p}_{2}},\,{{p}_{3}})$ из (17) совпадает со значениями $W({{p}_{1}},\,{{p}_{2}},\,{{p}_{3}})$ из (16). Далее через $\,\Pi $ и $\,{{\Pi }_{2}}$ будем обозначать целочисленные параллелепипеды с числом дискретных аргументов по каждой оси ${{N}_{1}},{\kern 1pt} {{N}_{2}},{{N}_{3}}$ и $2{{N}_{1}},{\kern 1pt} 2{{N}_{2}},2{{N}_{3}}$ соответственно. Теперь, проводя дискретное преобразование Фурье по каждой переменной от обеих частей (17), получим следующее равенство: Число арифметических операций и объем требуемой памяти для вычисления $W(p),\,\,p \in \Pi $, оцениваются формулами [5, 17]: В (19) ${\text{LOG}}(N){\kern 1pt} $ обозначает целочисленный логарифм, т.е. сумму всех простых делителей числа N. Например, ${\text{LOG}}(1000) = 21{\kern 1pt} $.

ЗАКЛЮЧЕНИЕ

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

Авторы заявляют об отсутствии конфликта интересов.