Радиотехника и электроника, 2019, T. 64, № 10, стр. 943-952

Две методики решения задачи дифракции плоской волны на теле вращения, расположенном в диэлектрическом полупространстве

С. А. Маненков *

Московский технический университет связи и информатики
111024 Москва, ул. Авиамоторная, 8а, Российская Федерация

* E-mail: mail44471@mail.ru

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

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

Аннотация

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

ВВЕДЕНИЕ

Задача дифракции волн на компактных препятствиях, расположенных в неоднородных средах представляет большой практический интерес. В частности, в литературе рассматривались различные двумерные задачи дифракции электромагнитных волн на бесконечных цилиндрических телах [14], а также задачи дифракции на телах вращения [57], расположенных в плоскослоистых средах.

В данной работе исследована задача рассеяния плоской волны на диэлектрическом теле вращения, погруженном в однородное диэлектрическое полупространство. Предполагается, что ось вращения тела перпендикулярна границе раздела сред. Для решения задачи применялись две различные методики. Первый метод основан на использовании модифицированного метода дискретных источников (ММДИ), с помощью которого ранее были рассмотрены задачи дифракции на одиночном теле вращения [8, 9], дифракции на группе тел [10, 11], рассеяние моды круглого диэлектрического волновода на теле вращения, расположенном внутри волновода [12], и др. При этом материальные характеристики среды тела предполагаются не зависящими от координат. Для решения задачи дифракции на рассеивателе, расположенном в полупространстве, применялась матричная функция Грина (ФГ), учитывающая слоистый характер среды. Такой подход позволяет свести задачу к системе интегральных уравнений (СИУ) относительно некоторых функций (токов), распределенных на вспомогательных поверхностях, расположенных внутри и вне исходной поверхности тела. При этом отпадает необходимость рассматривать токи на границе раздела сред.

Второй метод является приближенным. Назовем его модифицированным методом объемного интегрального уравнения (ММОИУ). В его основе лежит объемное интегральное уравнение (ИУ) 2-го рода относительно вектора напряженности электрического поля внутри тела. В отличие от первого метода здесь предполагается, что магнитная проницаемость во всех средах одинакова, а диэлектрическая проницаемость тела является функцией координат, не зависящей от угловой координаты (в цилиндрической системе координат, связанной с телом вращения). Отметим, что данная работа является обобщением работы [13] на случай дифракции на неоднородном теле вращения, погруженном в полупространство. В основе метода решения задачи лежат две идеи. Во-первых, как и в случае первой методики, использовалась матричная ФГ однородного полупространства. Это позволило свести задачу к объемному ИУ по области внутри тела. Во-вторых, ядро ИУ разбивали на сингулярную и регулярную (обусловленную границей раздела) части. Как и в работе [13] сингулярную часть ядра заменяли близким ядром, не имеющим особенности при совпадении аргументов. Алгебраизацию ИУ проводили при использовании метода Крылова–Боголюбова (с предварительным разложением неизвестной функции в ряд Фурье по угловой координате). В силу гладкости регулярной части ядра соответствующие интегралы, через которые выражаются добавочные матричные элементы алгебраической системы (обусловленные наличием регулярной части ядра), вычисляли при помощи формулы, являющейся двумерным аналогом формулы прямоугольников. Это позволило существенно сократить объем требуемых вычислений.

1. РЕШЕНИЕ ЗАДАЧИ ДИФРАКЦИИ НА ОДНОРОДНОМ ТЕЛЕ ВРАЩЕНИЯ ПРИ ПОМОЩИ ММДИ

Рассмотрим математическую постановку задачи. Пусть компактный рассеиватель в виде тела вращения расположен внутри диэлектрического полупространства с характеристиками ${{\varepsilon }_{2}},$ ${{\mu }_{2}}.$ Проницаемости среды второго полупространства (которое не содержит рассеивателя) обозначим через ${{\varepsilon }_{1}},$ ${{\mu }_{1}},$ характеристики среды внутри тела – ${{\varepsilon }_{i}},$ ${{\mu }_{i}}.$ Будем считать, что ось вращения тела перпендикулярна границе раздела сред. Введем декартову систему координат, причем ось $z$ направим вдоль оси вращения рассеивателя (рис. 1). Выберем начало системы координат внутри области $D$, занимаемой телом. Обозначим через $S$ границу области $D$. Пусть $d$ – расстояние от начала координат до границы раздела сред. Предположим, что на поверхности $S$ рассеивателя выполнены условия сопряжения:

(1)
${{\left. {\left[ {\vec {n} \times \vec {E}} \right]} \right|}_{S}} = 0,\,\,\,\,{{\left. {\left[ {\vec {n} \times \vec {H}} \right]} \right|}_{S}} = 0.$
Здесь $\vec {n}$ – нормаль, внешняя к поверхности тела, $\vec {E}$, $\vec {H}$ – полное электрическое и магнитное поле, квадратные скобки означают скачок соответствующей величины. На границе раздела сред $z = d$ выполнены следующие условия сопряжения:
(2)
${{\left. {\left[ {{{{\vec {i}}}_{z}} \times \vec {E}} \right]} \right|}_{{z = d}}} = 0,\,\,\,\,{{\left. {\left[ {{{{\vec {i}}}_{z}} \times \vec {H}} \right]} \right|}_{{z = d}}} = 0,$
где ${{\vec {i}}_{z}}$ – единичный орт вдоль оси z.

Рис. 1.

Геометрия задачи.

Предполагаем, что структура облучается плоской волной, падающей из верхнего полупространства ($z > d$), причем будем рассматривать две поляризации падающего поля. В первом случае

(3)
$E_{y}^{i} = \exp ( - i{{k}_{1}}x\sin {{\chi }_{1}} + i{{k}_{1}}z\cos {{\chi }_{1}}),\,\,\,\,E_{x}^{i} = E_{z}^{i} = 0,$
а во втором
(4)
$\begin{gathered} H_{y}^{i} = \exp ( - i{{k}_{1}}x\sin {{\chi }_{1}} + i{{k}_{1}}z\cos {{\chi }_{1}}), \\ H_{x}^{i} = H_{z}^{i} = 0, \\ \end{gathered} $
где ${{\chi }_{1}}$ – угол падения плоской волны, ${{k}_{1}} = \omega \sqrt {{{\varepsilon }_{1}}{{\mu }_{1}}} $ – волновое число среды верхнего полупространства, $\omega $ – круговая частота. Остальные компоненты падающего поля определяются из уравнений Максвелла. Введем также в рассмотрение первичное и вторичное (рассеянное) поля. А именно,
(5)
$\vec {E} = {{\vec {E}}^{0}} + {{\vec {E}}^{1}},\,\,\,\,\vec {H} = {{\vec {H}}^{0}} + {{\vec {H}}^{1}},$
где первичное поле ${{\vec {E}}^{0}},{{\vec {H}}^{0}}$ определяется из задачи дифракции плоской волны (3) или (4) на границе раздела сред в отсутствие компактного тела. Решение данной задачи общеизвестно и здесь не приводится. Будем также считать, что на бесконечности выполнены условия излучения для рассеянного поля:
(6)
$\mathop {\lim }\limits_{r \to \infty } {{\vec {E}}^{1}} = 0,\,\,\,\,\mathop {\lim }\limits_{r \to \infty } {{\vec {H}}^{1}} = 0,\,\,\,\,\operatorname{Im} {{\varepsilon }_{p}},\,\,\,{{\mu }_{p}} < 0,\,\,\,\,p = 1,2,$
где $r$ – расстояние от начала координат до точки наблюдения. Таким образом, в соответствии с принципом предельного поглощения вначале предполагаем, что проницаемости всех сред имеют малые мнимые части, а затем в окончательных формулах устремляем мнимые части к нулю.

В соответствии со схемой ММДИ, представим поле внутри и вне тела в виде:

(7)
$\vec {E}(\vec {r}) = \left\{ \begin{gathered} \int\limits_{{{\Sigma }_{2}}} {{\mathbf{G}}_{2}^{E}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{2}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} ,\,\,\,\,\vec {r} \in D, \hfill \\ {{{\vec {E}}}^{0}}(\vec {r}) + \int\limits_{{{\Sigma }_{1}}} {{\mathbf{G}}_{1}^{E}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{1}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} ,\,\,\,\,\vec {r} \in {{\mathbb{R}}^{3}}\,\backslash \,D, \hfill \\ \end{gathered} \right.$
(8)
$\vec {H}(\vec {r}) = \left\{ \begin{gathered} \int\limits_{{{\Sigma }_{2}}} {{\mathbf{G}}_{2}^{H}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{2}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} ,\,\,\,\vec {r} \in D, \hfill \\ {{{\vec {H}}}^{0}}(\vec {r}) + \int\limits_{{{\Sigma }_{1}}} {{\mathbf{G}}_{1}^{H}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{1}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} ,\,\,\,\,\vec {r} \in {{\mathbb{R}}^{3}}\,\backslash \,D, \hfill \\ \end{gathered} \right.$
где ${\mathbf{G}}_{1}^{{E,H}}$ и ${\mathbf{G}}_{2}^{{E,H}}$ – матричные ФГ слоистой среды вне рассеивателя и среды внутри тела, ${{\Sigma }_{1}}$и ${{\Sigma }_{2}}$ – вспомогательные поверхности вращения, расположенные внутри и вне исходной поверхности тела $S$, а ${{\vec {J}}_{1}}$ и ${{\vec {J}}_{2}}$ – неизвестные токи, распределенные на поверхностях ${{\Sigma }_{1}}$ и ${{\Sigma }_{2}}$.

Рассмотрим вопрос о выборе вспомогательных поверхностей ${{\Sigma }_{1}}$ и ${{\Sigma }_{2}}$. Предположим вначале, что поверхность S задана в сферических координатах:

(9)
$x = r\sin \theta \cos \varphi ,\,\,\,\,y = r\sin \theta \sin \varphi ,\,\,\,z = r\cos \theta ,$
где $r = r(\theta )$ и $\theta \in [0,\pi ].$ Тогда уравнения вспомогательных поверхностей в сферической системе координат имеют следующий вид [14]:
(10)
$\begin{gathered} {{r}_{{{{\Sigma }_{p}}}}} = \left| \xi \right|,\,\,\,\,{{\theta }_{{{{\Sigma }_{p}}}}} = \arg \xi , \\ \xi (\theta ) = r(\theta \pm i{{\delta }_{ \pm }})\exp (i\theta \mp {{\delta }_{ \pm }}),\,\,\,\,p = 1,2, \\ \end{gathered} $
где ${{r}_{{{{\Sigma }_{p}}}}},{{\theta }_{{{{\Sigma }_{p}}}}}$ – сферические координаты точки на вспомогательной поверхности ${{\Sigma }_{p}},$ ${{\delta }_{ \pm }}$ – положительные параметры, определяющие степень деформации исходной поверхности тела (знак плюс в (10) соответствует внутренней поверхности ${{\Sigma }_{1}}$, минус – внешней поверхности ${{\Sigma }_{2}}$). Указанные параметры изменяются в диапазоне $0 < {{\delta }_{ \pm }} < \delta _{{\max }}^{ \pm },$ где выбор величин $\delta _{{\max }}^{ \pm },$ которые определяют максимальные значения параметров ${{\delta }_{ \pm }},$ подробно рассмотрен в работе [14]. Для нахождения декартовых координат точки на вспомогательных поверхностях используем формулы
(11)
${{x}_{{{{\Sigma }_{p}}}}} = \operatorname{Im} \xi \cos \varphi ,\,\,\,\,{{y}_{{{{\Sigma }_{p}}}}} = \operatorname{Im} \xi \sin \varphi ,\,\,\,\,{{z}_{{{{\Sigma }_{p}}}}} = \operatorname{Re} \xi ,$
причем $\xi = {{z}_{{{{\Sigma }_{p}}}}} + i{{\rho }_{{{{\Sigma }_{p}}}}},$ где ${{\rho }_{{{{\Sigma }_{p}}}}}$ и ${{z}_{{{{\Sigma }_{p}}}}}$ – цилиндрические координаты точки на вспомогательной поверхности, $p = 1,2.$

Пусть далее поверхность тела задана в вытянутых сфероидальных координатах:

(12)
$\begin{gathered} x = f\operatorname{sh} \alpha \sin \beta \cos \varphi ,\,\,\,\,y = f\operatorname{sh} \alpha \sin \beta \sin \varphi , \\ z = f\operatorname{ch} \alpha \cos \beta , \\ \end{gathered} $
причем уравнение поверхности S имеет вид $\alpha = \alpha (\beta ),$ где $\beta \in [0,\pi ].$ Тогда вспомогательная поверхность определяется соотношениями [15]
(13)
$\begin{gathered} {{\alpha }_{{{{\Sigma }_{p}}}}} = \operatorname{Re} \eta ,\,\,\,\,{{\beta }_{{{{\Sigma }_{p}}}}} = \operatorname{Im} \eta , \\ \eta (\beta ) = \alpha (\beta \pm i{{\delta }_{ \pm }}) + i(\beta \pm i{{\delta }_{ \pm }}), \\ \end{gathered} $
где $\left( {{{\alpha }_{{{{\Sigma }_{p}}}}},{{\beta }_{{{{\Sigma }_{p}}}}},\varphi } \right)$ – сфероидальные координаты “образа” точки с координатами $\left( {\alpha ,\beta ,\varphi } \right),$ расположенной на поверхности S. Для получения декартовых координат точки на вспомогательной поверхности нужно вновь использовать формулы (11), в которых в данном случае $\xi (\beta ) = f\,{\text{ch}}\eta \,(\beta ).$ В случае сплюснутых сфероидальных координат $\xi (\beta ) = f\,{\text{sh}}\eta \,(\beta )$. Заметим, что в случае сферических координат мы обозначим $\alpha = \ln r,\,\,\beta = \theta .$ Таким образом, в дальнейшем будем считать, что поверхность тела задана уравнением $\alpha = \alpha (\beta )$ в соответствующей ортогональной системе координат.

Матричные ФГ ${\mathbf{G}}_{2}^{{E,H}}$ внутри тела определяются в соответствии с формулами [16]

(14)
${\mathbf{G}}_{2}^{E}(\vec {r},\vec {r}{\kern 1pt} '){{\vec {J}}_{2}}(\vec {r}{\kern 1pt} ') = {{\vec {A}}_{i}} + \frac{1}{{k_{i}^{2}}}\nabla \left( {\nabla \cdot {{{\vec {A}}}_{i}}} \right),$
(15)
${\mathbf{G}}_{2}^{H}(\vec {r},\vec {r}{\kern 1pt} '){{\vec {J}}_{2}}(\vec {r}{\kern 1pt} ') = \frac{i}{{{{k}_{i}}{{\zeta }_{i}}}}\nabla \times {{\vec {A}}_{i}},$
где ${{k}_{i}} = \omega \sqrt {{{\varepsilon }_{i}}{{\mu }_{i}}} $ – волновое число среды внутри тела, ${{\zeta }_{i}} = \sqrt {{{{{\mu }_{i}}} \mathord{\left/ {\vphantom {{{{\mu }_{i}}} {{{\varepsilon }_{i}}}}} \right. \kern-0em} {{{\varepsilon }_{i}}}}} $ – волновой импеданс среды тела,

(16)
$\begin{gathered} {{{\vec {A}}}_{i}} = {{{\vec {J}}}_{2}}(\vec {r}{\kern 1pt} '){{G}_{i}}(\vec {r},\vec {r}{\kern 1pt} '),\,\,\,\,{{G}_{i}} = \frac{{\exp ( - i{{k}_{i}}R)}}{{4\pi R}}, \\ R = \left| {\vec {r} - \vec {r}{\kern 1pt} '} \right|. \\ \end{gathered} $

Для построения ФГ слоистой среды, в которой расположено тело, т.е. матричных функций ${\mathbf{G}}_{1}^{E}$ и ${\mathbf{G}}_{1}^{H},$ применим методику, описанную в работе [17]. Для вывода представлений этих величин рассматривалась вспомогательная задача. Введем обозначения:

(17)
${{\vec {E}}^{G}}(\vec {r},\vec {r}{\kern 1pt} ') = {\mathbf{G}}_{1}^{E}(\vec {r},\vec {r}{\kern 1pt} '){{\vec {J}}_{1}}(\vec {r}{\kern 1pt} '),$
(18)
${{\vec {H}}^{G}}(\vec {r},\vec {r}{\kern 1pt} ') = {\mathbf{G}}_{1}^{H}(\vec {r},\vec {r}{\kern 1pt} '){{\vec {J}}_{1}}(\vec {r}{\kern 1pt} ').$

Поля $\vec {E}_{{}}^{G}$ и $\vec {H}_{{}}^{G}$ удовлетворяют уравнениям Максвелла

(19)
$\begin{gathered} \nabla \times {{{\vec {E}}}^{G}} = - i{{k}_{2}}{{\zeta }_{2}}{{{\vec {H}}}^{G}},\,\,\,\,\nabla \times {{{\vec {H}}}^{G}} = \frac{{i{{k}_{2}}}}{{{{\zeta }_{2}}}}{{{\vec {E}}}^{G}} + \\ + \,\,\frac{i}{{{{k}_{2}}{{\zeta }_{2}}}}{{{\vec {J}}}_{1}}(\vec {r}{\kern 1pt} ')\delta (\vec {r} - \vec {r}{\kern 1pt} '),\,\,\,\,z < d, \\ \end{gathered} $
(20)
$\nabla \times {{\vec {E}}^{G}} = - i{{k}_{1}}{{\zeta }_{1}}{{\vec {H}}^{G}},\,\,\,\,\nabla \times {{\vec {H}}^{G}} = \frac{{i{{k}_{1}}}}{{{{\zeta }_{1}}}}{{\vec {E}}^{G}},\,\,\,\,z > d,$
условиям сопряжения вида (2) на границе раздела и условиям (6) на бесконечности. В (19) и (20) ${{k}_{2}} = \omega \sqrt {{{\varepsilon }_{2}}{{\mu }_{2}}} $ – волновое число среды нижнего полупространства, ${{\zeta }_{1}} = \sqrt {{{{{\mu }_{1}}} \mathord{\left/ {\vphantom {{{{\mu }_{1}}} {{{\varepsilon }_{1}}}}} \right. \kern-0em} {{{\varepsilon }_{1}}}}} $ и ${{\zeta }_{2}} = \sqrt {{{{{\mu }_{2}}} \mathord{\left/ {\vphantom {{{{\mu }_{2}}} {{{\varepsilon }_{2}}}}} \right. \kern-0em} {{{\varepsilon }_{2}}}}} $ – волновые импедансы сред верхнего и нижнего полупространств. В дальнейшем перейдем к цилиндрическим координатам $(\rho ,\varphi ,z).$ Тогда для нахождения компонент матричных ФГ ${\mathbf{G}}_{1}^{E}$ и ${\mathbf{G}}_{1}^{H}$ необходимо рассмотреть три задачи, соответствующие возбуждению слоистой среды сторонними токами следующего вида:
(21)
$\begin{gathered} {\text{I}}.\,\,\,\,{{{\vec {J}}}_{1}} = {{{\vec {i}}}_{{\rho {\kern 1pt} '}}}, \hfill \\ {\text{II}}.\,\,\,{{{\vec {J}}}_{1}} = {{{\vec {i}}}_{{\varphi {\kern 1pt} '}}}, \hfill \\ {\text{III}}.\,\,{{{\vec {J}}}_{1}} = {{{\vec {i}}}_{{z{\kern 1pt} '}}}, \hfill \\ \end{gathered} $
где $(\rho {\kern 1pt} ',\varphi {\kern 1pt} ',z{\kern 1pt} ')$ – цилиндрические координаты точки, в которой расположен точечный источник, а ${{\vec {i}}_{{\rho {\kern 1pt} '}}},{{\vec {i}}_{{\varphi {\kern 1pt} '}}},{{\vec {i}}_{{z{\kern 1pt} '}}}$ – единичные орты по осям координат в точке источника. При этом координаты векторов $\vec {E}_{{}}^{G}$ и ${{\vec {H}}^{G}},$ полученные при решении задач I–III образуют столбцы матричных ФГ ${\mathbf{G}}_{1}^{E}$ и ${\mathbf{G}}_{1}^{H}.$

Во всех трех случаях удобно представить функции $\vec {E}_{{}}^{G}$ и $\vec {H}_{{}}^{G}$ в виде

(22)
${{\vec {E}}^{G}} = \vec {A} + \frac{1}{{{{k}^{2}}(z)}}\nabla (\nabla \cdot \vec {A}),$
(23)
${{\vec {H}}^{G}} = \frac{i}{{k(z)\zeta (z)}}\nabla \times \vec {A},$
где

(24)
$k(z) = \left\{ \begin{gathered} {{k}_{1}},\,\,z > d, \hfill \\ {{k}_{2}},\,\,z < d; \hfill \\ \end{gathered} \right.\,\,\,\,\zeta (z) = \left\{ \begin{gathered} {{\zeta }_{1}},\,\,z > d, \hfill \\ {{\zeta }_{2}},\,\,z < d. \hfill \\ \end{gathered} \right.$

В формулах (22) и (23) потенциал $\vec {A}$ определяется из соотношения (в декартовых координатах) [17]:

(25)
$\vec {A} = \left( {\begin{array}{*{20}{c}} {{{f}_{{11}}}}&0&0 \\ 0&{{{f}_{{11}}}}&0 \\ {\frac{{\partial {{f}_{{31}}}}}{{\partial x}}}&{\frac{{\partial {{f}_{{31}}}}}{{\partial y}}}&{{{f}_{{33}}}} \end{array}} \right){{\vec {J}}_{1}}.$

Функции ${{f}_{{11}}},{{f}_{{31}}},{{f}_{{33}}}$ представляются в виде

(26)
${{f}_{{11}}} = \frac{1}{{2\pi }}\int\limits_0^\infty {{{{v}}_{{11}}}(\kappa ,z,z{\kern 1pt} '){{J}_{0}}(\kappa {{R}_{ \bot }})\kappa d\kappa } ,$
(27)
${{f}_{{31}}} = \frac{1}{{2\pi }}\int\limits_0^\infty {{{v}_{{31}}}(\kappa ,z,z{\kern 1pt} '){{J}_{0}}(\kappa {{R}_{ \bot }})\kappa d\kappa } ,$
(28)
${{f}_{{33}}} = \frac{1}{{2\pi }}\int\limits_0^\infty {{{v}_{{33}}}(\kappa ,z,z{\kern 1pt} '){{J}_{0}}(\kappa {{R}_{ \bot }})\kappa d\kappa } ,$
где ${{R}_{ \bot }} = \sqrt {{{{(x - x{\kern 1pt} ')}}^{2}} + {{{(y - y{\kern 1pt} ')}}^{2}}} ,$ причем спектральные функции ${{v}_{{11}}},{{v}_{{31}}},{{v}_{{33}}}$ являются решениями одномерных краевых задач [17]

(29)
(30)
(31)

В формулах (29)(31) функции $\varepsilon (z)$ и $\mu (z)$ определяются аналогично соотношениям (24), штрих означает производную по переменной z. Решения задач (29)–(31) находятся аналогично тому, как это описано в работе [17]. Далее, нетрудно показать, что в цилиндрических координатах потенциал $\vec {A}$ выражается следующим образом:

(32)
$\begin{gathered} {{A}_{\rho }} = {{f}_{{11}}}\cos (\varphi - \varphi '),\,\,\,\,{{A}_{\varphi }} = - {{f}_{{11}}}\sin (\varphi - \varphi {\kern 1pt} '), \\ {{A}_{z}} = - \frac{{\partial {{f}_{{31}}}}}{{\partial \rho {\kern 1pt} '}}, \\ \end{gathered} $
(33)
$\begin{gathered} {{A}_{\rho }} = {{f}_{{11}}}\sin (\varphi - \varphi {\kern 1pt} '),\,\,\,\,{{A}_{\varphi }} = {{f}_{{11}}}\cos (\varphi - \varphi {\kern 1pt} '), \\ {{A}_{z}} = - \frac{1}{{\rho {\kern 1pt} '}}\frac{{\partial {{f}_{{31}}}}}{{\partial \varphi {\kern 1pt} '}}, \\ \end{gathered} $
(34)
${{A}_{\rho }} = 0,\,\,\,{{A}_{\varphi }} = 0,\,\,\,{{A}_{z}} = {{f}_{{33}}},$
соответственно для задач I, II и III.

Перейдем к решению поставленной задачи дифракции. Для этого подставим формулы (7) и (8) в граничные условия (1). В результате получим следующую СИУ

(35)
$\begin{gathered} \vec {n} \times \int\limits_{{{\Sigma }_{1}}} {{\mathbf{G}}_{1}^{E}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{1}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} - \\ - \,\,\vec {n} \times \int\limits_{{{\Sigma }_{2}}} {{\mathbf{G}}_{2}^{E}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{2}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} = - \vec {n} \times {{{\vec {E}}}^{0}}, \\ \end{gathered} $
(36)
$\begin{gathered} \vec {n} \times \int\limits_{{{\Sigma }_{1}}} {{\mathbf{G}}_{1}^{H}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{1}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} - \vec {n} \times \int\limits_{{{\Sigma }_{2}}} {{\mathbf{G}}_{2}^{H}(\vec {r},\vec {r}{\kern 1pt} '){{{\vec {J}}}_{2}}(\vec {r}{\kern 1pt} ')ds{\kern 1pt} '} = \\ = - \vec {n} \times {{{\vec {H}}}^{0}},\,\,\,\,\vec {r} \in S. \\ \end{gathered} $

Пусть, как указано выше, поверхность тела S задана уравнением $\alpha = \alpha (\beta ),$ а параметрические уравнения вспомогательных поверхностей ${{\Sigma }_{1}}$ и ${{\Sigma }_{2}}$ имеют вид $\alpha = {{\alpha }_{{{{\Sigma }_{p}}}}}(\beta {\kern 1pt} '),$ $\beta = {{\beta }_{{{{\Sigma }_{p}}}}}(\beta {\kern 1pt} ')$ ($p = 1,2$), причем $\beta {\kern 1pt} ' \in [0,\pi ]$ (см. формулы (13)). Для алгебраизации задачи разложим неизвестные функции (токи) и матричные ФГ в ряды Фурье:

(37)
${{\vec {J}}_{p}}({{\alpha }_{{{{\Sigma }_{p}}}}},{{\beta }_{{{{\Sigma }_{p}}}}},\varphi ) = \sum\limits_{m = - \infty }^\infty {{{{\vec {J}}}_{{p,m}}}({{\alpha }_{{{{\Sigma }_{p}}}}},{{\beta }_{{{{\Sigma }_{p}}}}})\exp (im\varphi )} ,$
(38)
$\begin{gathered} {\mathbf{G}}_{p}^{E}(\alpha ,\beta ,{{\alpha }_{{{{\Sigma }_{p}}}}},{{\beta }_{{{{\Sigma }_{p}}}}},\psi ) = \\ = \sum\limits_{m = - \infty }^\infty {{\mathbf{G}}_{{p,m}}^{E}(\alpha ,\beta ,{{\alpha }_{{{{\Sigma }_{p}}}}},{{\beta }_{{{{\Sigma }_{p}}}}})\exp (im\psi )} , \\ \end{gathered} $
(39)
$\begin{gathered} {\mathbf{G}}_{p}^{H}(\alpha ,\beta ,{{\alpha }_{{{{\Sigma }_{p}}}}},{{\beta }_{{{{\Sigma }_{p}}}}},\psi ) = \\ = \sum\limits_{m = - \infty }^\infty {{\mathbf{G}}_{{p,m}}^{H}(\alpha ,\beta ,{{\alpha }_{{{{\Sigma }_{p}}}}},{{\beta }_{{{{\Sigma }_{p}}}}})\exp (im\psi )} , \\ \psi = \varphi - \varphi {\text{'}},\,\,\,\,p = 1,2. \\ \end{gathered} $

Сделаем замену неизвестных гармоник тока [15]:

(40)
${{\vec {J}}_{{p,m}}} = {{\vec {I}_{p}^{m}} \mathord{\left/ {\vphantom {{\vec {I}_{p}^{m}} {\left( {{{h}_{\varphi }}{{h}_{{{{\Sigma }_{p}}}}}\sqrt {\dot {\alpha }_{{{{\Sigma }_{p}}}}^{2} + \dot {\beta }_{{{{\Sigma }_{p}}}}^{2}} } \right)}}} \right. \kern-0em} {\left( {{{h}_{\varphi }}{{h}_{{{{\Sigma }_{p}}}}}\sqrt {\dot {\alpha }_{{{{\Sigma }_{p}}}}^{2} + \dot {\beta }_{{{{\Sigma }_{p}}}}^{2}} } \right)}},$
где

(41)
$\vec {I}_{p}^{m} = - {{n}_{{{{\beta }_{p}}}}}I_{p}^{{m1}}{{\vec {i}}_{{{{\alpha }_{p}}}}} + {{n}_{{{{\alpha }_{p}}}}}I_{p}^{{m1}}{{\vec {i}}_{{{{\beta }_{p}}}}} + I_{p}^{{m2}}{{\vec {i}}_{\varphi }},\,\,\,\,p = 1,2.$

В формулах (40) и (41) ${{h}_{{{{\Sigma }_{p}}}}} = {{h}_{{{{\alpha }_{p}}}}} = {{h}_{{{{\beta }_{p}}}}}$ и ${{h}_{\varphi }}$ – коэффициенты Ламе, ${{\vec {i}}_{{{{\alpha }_{p}}}}},{{\vec {i}}_{{{{\beta }_{p}}}}},{{\vec {i}}_{\varphi }}$ – единичные орты выбранной системы координат, ${{n}_{{{{\alpha }_{p}}}}},{{n}_{{{{\beta }_{p}}}}}$ – координаты нормалей к поверхности ${{\Sigma }_{p}}$:

(42)
${{n}_{{{{\alpha }_{p}}}}} = \frac{{{{{\dot {\beta }}}_{{{{\Sigma }_{p}}}}}}}{{\sqrt {\dot {\alpha }_{{{{\Sigma }_{p}}}}^{2} + \dot {\beta }_{{{{\Sigma }_{p}}}}^{2}} }},\,\,\,\,{{n}_{{{{\beta }_{p}}}}} = - \frac{{{{{\dot {\alpha }}}_{{{{\Sigma }_{p}}}}}}}{{\sqrt {\dot {\alpha }_{{{{\Sigma }_{p}}}}^{2} + \dot {\beta }_{{{{\Sigma }_{p}}}}^{2}} }}.$

В формулах (40) и (42) через ${{\dot {\alpha }}_{{{{\Sigma }_{p}}}}}$ и ${{\dot {\beta }}_{{{{\Sigma }_{p}}}}}$ обозначены следующие производные:

(43)
$\begin{gathered} {{{\dot {\alpha }}}_{{{{\Sigma }_{p}}}}}(\beta {\kern 1pt} ') = \operatorname{Re} (\dot {\alpha }(\beta {\kern 1pt} '\,\, \pm i{{\delta }_{ \pm }}) + i), \\ {{{\dot {\beta }}}_{{{{\Sigma }_{p}}}}}(\beta {\kern 1pt} ') = \operatorname{Im} (\dot {\alpha }(\beta {\kern 1pt} '\,\, \pm i{{\delta }_{ \pm }}) + i), \\ \end{gathered} $

где точка в правых частях равенств (43) означает производную функции $\alpha (\beta {\kern 1pt} '\, \pm i{{\delta }_{ \pm }})$ по аргументу $\beta {\kern 1pt} '\, \pm i{{\delta }_{ \pm }},$ а в левых частях точка – производная по параметру $\beta {\kern 1pt} '.$ Верхний знак в соотношениях (43) соответствует внутренней поверхности ${{\Sigma }_{1}}.$

Далее подставим формулы (37)(39) в систему уравнений (35) и (36) и спроектируем полученные равенства на базис Фурье по угловой координате $\varphi $. Принимая во внимание соотношения (40) и (41), получим следующую бесконечную систему одномерных интегральных уравнений:

(44)
$\begin{gathered} {\mathbf{\hat {K}}}_{{11}}^{m}{\mathbf{I}}_{1}^{m} + {\mathbf{\hat {K}}}_{{12}}^{m}{\mathbf{I}}_{2}^{m} = {\mathbf{B}}_{1}^{m}, \\ {\mathbf{\hat {K}}}_{{21}}^{m}{\mathbf{I}}_{1}^{m} + {\mathbf{\hat {K}}}_{{22}}^{m}{\mathbf{I}}_{2}^{m} = {\mathbf{B}}_{2}^{m},\,\,\,\,m = 0, \pm 1, \pm 2,..., \\ \end{gathered} $
где ${\mathbf{I}}_{p}^{m} = {{\left[ {I_{p}^{{m1}},I_{p}^{{m2}}} \right]}^{T}},$ ${\mathbf{B}}_{p}^{m} = {{\left[ {B_{p}^{{m1}},B_{p}^{{m2}}} \right]}^{T}},$ $p = 1,2.$ В системе (44) вводятся следующие матричные операторы
(45)
$\begin{gathered} {\mathbf{\hat {K}}}_{{pq}}^{m}{\mathbf{I}}_{q}^{m} = \left( {\begin{array}{*{20}{c}} {\int\limits_0^\pi {K_{{pq,11}}^{m}(\beta ,\beta {\kern 1pt} ')I_{q}^{{m1}}(\beta {\kern 1pt} ')d\beta {\kern 1pt} '} }&{\int\limits_0^\pi {K_{{pq,12}}^{m}(\beta ,\beta {\kern 1pt} ')I_{q}^{{m2}}(\beta {\kern 1pt} ')d\beta {\kern 1pt} '} } \\ {\int\limits_0^\pi {K_{{pq,21}}^{m}(\beta ,\beta {\kern 1pt} ')I_{q}^{{m1}}(\beta {\kern 1pt} ')d\beta {\kern 1pt} '} }&{\int\limits_0^\pi {K_{{pq,22}}^{m}(\beta ,\beta {\kern 1pt} ')I_{q}^{{m2}}(\beta {\kern 1pt} ')d\beta {\kern 1pt} '} } \end{array}} \right), \\ p,q = 1,2. \\ \end{gathered} $
Ядра в (45) имеют вид
(46)
$\begin{gathered} K_{{11,11}}^{m} = {{n}_{\alpha }}({{n}_{{{{\alpha }_{1}}}}}e_{{1,22}}^{m} - {{n}_{{{{\beta }_{1}}}}}e_{{1,21}}^{m}) - {{n}_{\beta }}({{n}_{{{{\alpha }_{1}}}}}e_{{1,12}}^{m} - {{n}_{{{{\beta }_{1}}}}}e_{{1,11}}^{m}), \\ K_{{11,12}}^{m} = {{n}_{\alpha }}e_{{1,23}}^{m} - {{n}_{\beta }}e_{{1,13}}^{m},\,\,\,\,K_{{11,21}}^{m} = {{n}_{{{{\alpha }_{1}}}}}e_{{1,32}}^{m} - {{n}_{{{{\beta }_{1}}}}}e_{{1,31}}^{m}, \\ K_{{11,22}}^{m} = e_{{1,33}}^{m}, \\ \end{gathered} $
(47)
$\begin{gathered} K_{{12,11}}^{m} = - ({{n}_{\alpha }}({{n}_{{{{\alpha }_{2}}}}}e_{{2,22}}^{m} - {{n}_{{{{\beta }_{2}}}}}e_{{2,21}}^{m}) - {{n}_{\beta }}({{n}_{{{{\alpha }_{2}}}}}e_{{2,12}}^{m} - {{n}_{{{{\beta }_{2}}}}}e_{{2,11}}^{m})), \\ K_{{12,12}}^{m} = - ({{n}_{\alpha }}e_{{2,23}}^{m} - {{n}_{\beta }}e_{{2,13}}^{m}), \\ K_{{12,21}}^{m} = - ({{n}_{{{{\alpha }_{2}}}}}e_{{2,32}}^{m} - {{n}_{{{{\beta }_{2}}}}}e_{{2,31}}^{m}),\,\,K_{{12,22}}^{m} = - e_{{2,33}}^{m}, \\ \end{gathered} $
где ${{n}_{\alpha }},\,\,{{n}_{\beta }}$ – координаты нормали к поверхности S. Формулы для остальных ядер в СИУ (44) получаются заменой $e_{{p,jl}}^{m}$ на $h_{{p,jl}}^{m}$ в равенствах (46), (47). Величины $e_{{p,jl}}^{m}$ и $h_{{p,jl}}^{m}$ представляют собой элементы матричной функции ${\mathbf{G}}_{{p,m}}^{E}$ или ${\mathbf{G}}_{{p,m}}^{H},$ заданной в соответствующих ортогональных координатах. Эти величины наиболее просто могут быть выражены в цилиндрических координатах. При этом удобно представить ${\mathbf{G}}_{{1,m}}^{E}$ и ${\mathbf{G}}_{{1,m}}^{H}$ в виде суммы:
(48)
$\begin{gathered} {\mathbf{G}}_{{1,m}}^{E} = {\mathbf{G}}_{{1,m}}^{{E,s}} + {\mathbf{G}}_{{1,m}}^{{E,r}}, \\ {\mathbf{G}}_{{1,m}}^{H} = {\mathbf{G}}_{{1,m}}^{{H,s}} + {\mathbf{G}}_{{1,m}}^{{H,r}}, \\ \end{gathered} $
где ${\mathbf{G}}_{{1,m}}^{{E,s}},\;{\mathbf{G}}_{{1,m}}^{{H,s}}$ – части ФГ, соответствующие полю точечного источника в однородной среде с волновым числом ${{k}_{2}}$ (сингулярные части ФГ, определяемые аналогично формулам (14)(16)), и ${\mathbf{G}}_{{1,m}}^{{E,r}},\;{\mathbf{G}}_{{1,m}}^{{H,r}}$ – добавочные части ФГ, обусловленные наличием границы раздела сред. Подробные выражения для элементов ФГ ${\mathbf{G}}_{{1,m}}^{{E,s}},\,\,{\mathbf{G}}_{{1,m}}^{{H,s}},$ а также функций ${\mathbf{G}}_{{2,m}}^{E},\,\,{\mathbf{G}}_{{2,m}}^{H}$ приведены в работе [15]. Элементы добавочных ФГ ${\mathbf{G}}_{{1,m}}^{{E,r}},\,\,{\mathbf{G}}_{{1,m}}^{{H,r}}$ получаются из формул (17)–(34).

Для перехода от цилиндрических координат к ортогональным координатам $(\alpha ,\beta ,\varphi )$ используем соотношение [15]

(49)
${\mathbf{G}} = {{{\mathbf{{\rm T}}}}^{T}}{{{\mathbf{G}}}_{{{\text{cyl}}}}}{\mathbf{T}}{\kern 1pt} ',$
где ${{{\mathbf{G}}}_{{{\text{cyl}}}}}$ – тензор ${\mathbf{G}}_{{p,m}}^{E}$ или ${\mathbf{G}}_{{p,m}}^{H},$ заданный в цилиндрических координатах, а ${\mathbf{G}}$ – тот же тензор, заданный в указанных выше ортогональных координатах, ${\mathbf{{\rm T}}},\,\,{\mathbf{T}}{\kern 1pt} '$ – матрицы перехода. Матрица ${\mathbf{{\rm T}}}$ определяется следующим образом:
(50)
${\mathbf{T}} = \left( {\begin{array}{*{20}{c}} {\frac{1}{{\left| {\xi {\kern 1pt} '} \right|}}\operatorname{Im} \xi {\kern 1pt} '}&{\frac{1}{{\left| {\xi {\kern 1pt} '} \right|}}\operatorname{Re} \xi {\kern 1pt} '}&0 \\ 0&0&1 \\ {\frac{1}{{\left| {\xi {\kern 1pt} '} \right|}}\operatorname{Re} \xi {\kern 1pt} '}&{ - \frac{1}{{\left| {\xi {\kern 1pt} '} \right|}}\operatorname{Im} \xi {\kern 1pt} '}&0 \end{array}} \right),$
где штрих означает производную от $\xi $ по переменой $\eta = \alpha + i\beta .$ Матрица ${\mathbf{T}}{\kern 1pt} '$ получается из (50) путем замены$(\alpha ,\beta )$ на $({{\alpha }_{{{{\Sigma }_{1}}}}},{{\beta }_{{{{\Sigma }_{1}}}}})$ или $({{\alpha }_{{{{\Sigma }_{2}}}}},{{\beta }_{{{{\Sigma }_{2}}}}}).$

2. РЕШЕНИЕ ЗАДАЧИ ДИФРАКЦИИ НА ТЕЛЕ ВРАЩЕНИЯ С ПЕРЕМЕННЫМИ МАТЕРИАЛЬНЫМИ ХАРАКТЕРИСТИКАМИ ПРИ ПОМОЩИ ММОИУ

Рассмотрим задачу дифракции поля плоской волны вида (3) или (4) на теле вращения, расположенном в однородном полупространстве. Геометрия задачи вновь иллюстрируется рис. 1. Как и выше, будем считать, что диэлектрические проницаемости верхнего и нижнего полупространств равны ${{\varepsilon }_{1}}$ и ${{\varepsilon }_{2}},$ а ${{\varepsilon }_{i}} = {{\varepsilon }_{i}}(\rho ,z)$ – диэлектрическая проницаемость внутри тела, которая не зависит от угла $\varphi $ (система координат та же, что и в разд. 1). В отличие от задачи, рассмотренной в разд. 1, будем предполагать среду во всем пространстве магнитооднородной, причем всюду магнитная проницаемость $\mu = 1.$ Граничные условия на поверхности тела $S$, на границе раздела $z = d$ и условия излучения остаются прежними.

Используя векторный аналог формулы Грина [16], нетрудно показать, что вектор напряженности рассеянного электрического поля магнитооднородной слоистой среды представляется в виде

(51)
${{\vec {E}}^{1}}(\vec {r}) = \int\limits_D {(k_{i}^{2}(\vec {r}{\kern 1pt} ') - k_{2}^{2}){\mathbf{G}}_{1}^{E}(\vec {r},\vec {r}{\kern 1pt} ')\vec {E}(\vec {r}{\kern 1pt} ')dV{\kern 1pt} '} .$

Поместим далее точку наблюдения в область $D$, в результате чего получим, что полное электрическое поле удовлетворяет следующему ИУ 2-го рода [16]:

(52)
$\begin{gathered} \vec {E}(\vec {r}) + \int\limits_D {(k_{2}^{2} - k_{i}^{2}(\vec {r}{\kern 1pt} ')){\mathbf{G}}_{1}^{E}(\vec {r},\vec {r}{\kern 1pt} ')\vec {E}(\vec {r}{\kern 1pt} ')dV{\kern 1pt} '} = {{{\vec {E}}}^{0}}(\vec {r}), \\ \vec {r} \in D, \\ \end{gathered} $
где ${{k}_{i}} = \omega \sqrt {{{\varepsilon }_{i}}(\vec {r}{\kern 1pt} ')} $ – волновое число среды внутри тела. При этом ФГ ${\mathbf{G}}_{1}^{E}$ разбивается на сумму сингулярной ${\mathbf{G}}_{1}^{{E,s}}$ и регулярной ${\mathbf{G}}_{1}^{{E,r}}$ частей, причем ${\mathbf{G}}_{1}^{{E,s}}$ определяется через скалярную ФГ

(53)
$G(\vec {r},\vec {r}{\kern 1pt} ') = \frac{{\exp ( - i{{k}_{2}}R)}}{{4\pi R}}.$

Основная идея предлагаемого подхода, как и в работе [13], состоит в замене ФГ свободного пространства $G(\vec {r},\vec {r}{\kern 1pt} ')$ на функцию

(54)
${{G}_{\delta }}(\vec {r},\vec {r}{\kern 1pt} ') = \frac{{\exp ( - i{{k}_{2}}{{R}_{\delta }})}}{{4\pi {{R}_{\delta }}}},$
где ${{R}_{\delta }} = \sqrt {{{\rho }^{2}} + \rho {\kern 1pt} {{'}^{2}} - 2\rho \rho {\kern 1pt} '\cos \psi + {{{(z - z{\kern 1pt} ')}}^{2}} + {{\delta }^{2}}} ,$ $\delta > 0$ – малый параметр. Таким образом, ядро уравнения (52) не будет иметь особенностей при совпадении аргументов. Очевидно, что в силу малости параметра $\delta $ заменяем ядро ИУ (52) близкой функцией. В результате получаем приближенное решение исходной краевой задачи, которое будет более точным при меньшем значении параметра $\delta $.

Как и в предыдущем разделе, для алгебраизации задачи разложим неизвестную функцию, первичное поле и ФГ в ряды Фурье и подставим эти разложения в уравнение (52). В результате получим бесконечный набор двумерных ИУ для гармоник Фурье вектора напряженности электрического поля:

(55)
$\begin{gathered} {{{\vec {E}}}_{m}}(\vec {\rho }) + \\ + \,\,\int\limits_{\tilde {D}} {(k_{2}^{2} - k_{i}^{2}(\vec {\rho }{\kern 1pt} ')){\mathbf{G}}_{{1,m}}^{E}(\vec {\rho },\vec {\rho }{\kern 1pt} '){{{\vec {E}}}_{m}}(\vec {\rho }{\kern 1pt} ')\rho {\kern 1pt} 'd\rho {\kern 1pt} 'dz{\kern 1pt} '} = \\ = \vec {E}_{m}^{0}(\vec {\rho }),\,\,\,\,\vec {\rho } \in \tilde {D}. \\ \end{gathered} $

Здесь $m = 0, \pm 1, \pm 2...,$ $\vec {\rho } = (\rho ,z),$ $\tilde {D}$ – двумерная область, получаемая в осевом сечении области $D$ тела, ${{\vec {E}}_{m}}(\vec {\rho })$ – гармоники электрического поля. Для численного решения уравнений (55) будем предполагать, что область $\tilde {D}$ задана уравнениями

(56)
$\rho = \rho (u,v),\,\,\,\,z = z(u,v),$

где ${{u}_{1}} \leqslant u \leqslant {{u}_{2}},$ ${{v}_{1}} \leqslant v \leqslant {{v}_{2}}.$ Таким образом, имеется отображение области $\tilde {D}$ на прямоугольник $\tilde {D}{\kern 1pt} ' = [{{u}_{1}},{{u}_{2}}] \times [{{v}_{1}},{{v}_{2}}].$ В новых переменных ИУ (55) могут быть записаны в виде

(57)
$\begin{gathered} {{{\vec {E}}}_{m}}(u,v) + \int\limits_{{{u}_{1}}}^{{{u}_{2}}} {\int\limits_{{{v}_{1}}}^{{{{v}}_{2}}} {(k_{2}^{2} - k_{i}^{2}(u{\kern 1pt} ',v{\kern 1pt} ')){\mathbf{G}}_{{1,m}}^{E}(u,v,u{\kern 1pt} ',v{\kern 1pt} ')} } \times \\ \times \,\,{{{\vec {E}}}_{m}}(u{\kern 1pt} ',v{\kern 1pt} ')\left| {J(u{\kern 1pt} ',v{\kern 1pt} ')} \right|\rho {\kern 1pt} '(u{\kern 1pt} ',v{\kern 1pt} ')du{\kern 1pt} 'dv{\kern 1pt} ' = \\ = \vec {E}_{m}^{0}(u,v), \\ \end{gathered} $

где $m = 0, \pm 1, \pm 2...$, $J(u{\kern 1pt} ',v{\kern 1pt} ')$ – якобиан преобразования (56). Для решения уравнений (57) использовали метод Крылова–Боголюбова, аналогично работе [13]. В результате применения данного метода система уравнений (57) сводится к бесконечному набору линейных алгебраических систем относительно неизвестных коэффициентов $\vec {c}_{{jl}}^{m}$ разложения величин ${{\vec {E}}_{m}}(u,v)$ по базису, состоящему из кусочно-постоянных функций [13]:

(58)
$\begin{gathered} \sum\limits_{j = 1}^{{{N}_{1}}} {\sum\limits_{l = 1}^{{{N}_{2}}} {{\mathbf{G}}_{{st}}^{{jl,m}}\vec {c}_{{jl}}^{m}} } = \vec {w}_{{st}}^{m},\,\,\,\,s = 1,2,...,{{N}_{1}}, \\ t = 1,2,...,{{N}_{2}}. \\ \end{gathered} $

Здесь

(59)
$\begin{gathered} {\mathbf{G}}_{{st}}^{{jl,m}} = {{\delta }_{{sj}}}{{\delta }_{{tl}}} + \int\limits_{{{{\bar {u}}}_{j}} - {{\Delta }_{1}}}^{{{{\bar {u}}}_{j}} + {{\Delta }_{1}}} {\int\limits_{{{{\bar {v}}}_{l}} - {{\Delta }_{2}}}^{{{{\bar {v}}}_{l}} + {{\Delta }_{2}}} {(k_{2}^{2} - k_{i}^{2}(u,v))} } \times \\ \times \,\,{\mathbf{G}}_{{1,m}}^{E}({{{\bar {u}}}_{s}},{{{\bar {v}}}_{t}},u,v)\left| {J(u,v)} \right|\rho (u,v)dudv{\text{,}} \\ \end{gathered} $
(60)
$\vec {w}_{{st}}^{m} = \vec {E}_{m}^{0}({{\bar {u}}_{s}},{{\bar {v}}_{t}}),$
${{\delta }_{{mn}}}$ – символ Кронекера. В системах (58) точки коллокации определяются из соотношений
(61)
${{\bar {u}}_{s}} = {{u}_{1}} + (2s - 1){{\Delta }_{1}},\,\,\,\,s = 1,2,...,{{N}_{1}},$
(62)
${{\bar {v}}_{t}} = {{v}_{1}} + (2t - 1){{\Delta }_{2}},\,\,\,\,t = 1,2,...,{{N}_{2}},$
где ${{\Delta }_{1}} = \frac{{{{u}_{2}} - {{u}_{1}}}}{{2{{N}_{1}}}},$ ${{\Delta }_{2}} = \frac{{{{v}_{2}} - {{v}_{1}}}}{{2{{N}_{2}}}}.$

Матричные элементы (59) были вычислены с помощью следующего алгоритма. Воспользуемся тем, что ФГ ${\mathbf{G}}_{{1,m}}^{E}$ имеет вид (48). При этом величины ${\mathbf{G}}_{{1,m}}^{{E,r}}$ являются медленно меняющимися функциями координат, а сингулярная часть ${\mathbf{G}}_{{1,m}}^{{E,s}}$ меняется очень быстро при приближении точки наблюдения к точке источника. Тогда интегралы в (59) можно найти по приближенной формуле

(63)
$\begin{gathered} {\mathbf{G}}_{{st}}^{{jl,m}} = {{\delta }_{{sj}}}{{\delta }_{{tl}}} + \int\limits_{{{{\bar {u}}}_{j}} - {{\Delta }_{1}}}^{{{{\bar {u}}}_{j}} + {{\Delta }_{1}}} {\int\limits_{{{{\bar {v}}}_{l}} - {{\Delta }_{2}}}^{{{{\bar {v}}}_{l}} + {{\Delta }_{2}}} {(k_{2}^{2} - k_{i}^{2}(u,v))} } \times \\ \times \,\,{\mathbf{G}}_{{1,m}}^{{E,s}}({{{\bar {u}}}_{s}},{{{\bar {v}}}_{t}},u,v)\left| {J(u,v)} \right|\rho (u,v)dudv + \\ + \,\,4(k_{2}^{2} - k_{i}^{2}({{{\bar {u}}}_{j}},{{{\bar {v}}}_{l}})){\mathbf{G}}_{{1,m}}^{{E,r}}({{{\bar {u}}}_{s}},{{{\bar {v}}}_{t}},{{{\bar {u}}}_{j}},{{{\bar {v}}}_{l}}) \times \\ \times \,\,\left| {J({{{\bar {u}}}_{j}},{{{\bar {v}}}_{l}})} \right|\rho ({{{\bar {u}}}_{j}},{{{\bar {v}}}_{l}}){{\Delta }_{1}}{{\Delta }_{2}}. \\ \end{gathered} $

Как следует из формулы (63) и формул для ФГ разд. 1, добавочные матричные элементы, обусловленные регулярной частью ФГ, выражаются в виде однократных интегралов (интегралов Зоммерфельда). Кроме того, аналогично работе [13] повторные интегралы в (63) следует вычислять более точно только для диагональных элементов матрицы. Интегралы во внедиагональных элементах можно вычислять, например, применяя формулы Гаусса с небольшим числом узлов. При таком подходе удается существенно сократить время вычисления матричных элементов на ЭВМ.

3. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Рассмотрим результаты численного моделирования. На рис. 2а и 2б представлены зависимости невязки краевых условий (1) на контуре осевого сечения однородной сферы и однородного суперэллипсоида вращения, расположенных в диэлектрическом полупространстве. Было рассматрено падение волны (3), угол ${{\chi }_{1}} = 0.$ Задачу решали с помощью метода, основанного на ММДИ (первая методика). В случае решения задачи дифракции на суперэллипсоиде использовали сферические координаты для построения вспомогательных поверхностей. СИУ (44) решали методом коллокации [812, 15]. Уравнение контура осевого сечения суперэллипсоида имело вид

(64)
${{\left( {\frac{\rho }{a}} \right)}^{{2q}}} + {{\left( {\frac{z}{b}} \right)}^{{2q}}} = 1.$
Рис. 2.

Распределение невязки на контуре осевого сечения сферы (а) и суперэллипсоида вращения (б).

Параметры задачи были следующие: радиус сферы ${{k}_{0}}a = 1.5,$ полуоси суперэллипсоида ${{k}_{0}}a = 1.5$ и ${{k}_{0}}b = 1.5,$ параметр $q = 20,$ диэлектрические проницаемости сред – ${{\varepsilon }_{1}} = 1,$ ${{\varepsilon }_{2}} = 2,$ ${{\varepsilon }_{i}} = 4,$ расстояние ${{k}_{0}}d = 2,$ где ${{k}_{0}} = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } \lambda }} \right. \kern-0em} \lambda }$ ($\lambda $ – длина волны). Магнитная проницаемость всюду была равна единице. Угол наблюдения $\varphi = 0.$ На рис. 2а число дискретных источников равно $N = 30$ (кривая 1) и $N = 50$ (кривая 2), а на рис. 2б – $N = 250$ (кривая 1) и $N = 350$ (кривая 2). Как видно из рисунков, максимальный уровень невязки становится существенно меньше при увеличении числа дискретных источников. При этом максимальное значение модуля невязки не превосходит $5 \times {{10}^{{ - 10}}}$ в случае дифракции на сфере (при $N = 50$) и не превосходит $5 \times {{10}^{{ - 4}}}$ в случае дифракции на суперэллипсоиде (при $N = 350$). Данный факт гарантирует корректность методики, основанной на ММДИ.

В табл. 1 приведены значения модуля диаграммы рассеяния $\vec {F}(\theta ,\varphi )$ в верхнем полупространстве для рассматриваемой задачи дифракции, полученные при помощи обеих методик. Диаграмма рассеяния определяется по формуле

(65)
${{\vec {E}}^{1}} \approx \vec {F}(\theta ,\varphi )\frac{{\exp ( - i{{k}_{1}}r)}}{r},$
Таблица 1.  

Сравнение результатов расчета диаграммы рассеяния однородной сферы и однородного диэлектрического кругового цилиндра, полученных двумя методами

Угол, град Дифракция на сфере Дифракция на цилиндре
метод на основе ММДИ метод на основе ММОИУ метод на основе ММДИ метод на основе ММОИУ
0 0.2520 0.2542 0.5376 0.5462
9 0.2487 0.2508 0.5328 0.5413
18 0.2387 0.2408 0.5181 0.5264
27 0.2225 0.2245 0.4934 0.5013
36 0.2007 0.2024 0.4583 0.4656
45 0.1741 0.1756 0.4128 0.4193
54 0.1440 0.1453 0.3570 0.3626
63 0.1118 0.1128 0.2911 0.2955
72 0.07824 0.07896 0.2137 0.2169
81 0.04243 0.04282 0.1202 0.1220
90 0.1884 × 10–16 0.1901 × 10–16 0.5409 × 10–16 0.5490 × 10–16

где $r \to \infty .$ Конкретный вид диаграммы можно получить из нижней формулы в (7) в случае применения первой методики и из формулы (51) в случае использования второй методики (с помощью метода перевала). Отметим, что в последнем случае можно вычислять кратные интегралы, через которые выражается диаграмма рассеяния, применяя двумерный аналог формулы прямоугольников. Таким образом, нахождение диаграммы (после определения токовых функций) не требует большого объема вычислений, как в случае первой, так и в случае второй методики.

Результаты табл. 1 приведены для задачи рассеяния плоской волны (3) на однородной сфере и однородном цилиндре конечных размеров, погруженных в однородное полупространство. При расчете дифракции на цилиндре в случае использования метода на основе ММДИ поверхность цилиндра аппроксимировалась поверхностью суперэллипсоида (64) с параметром $q = 20.$ Радиус сферы ${{k}_{0}}a = 1.5,$ высота цилиндра ${{k}_{0}}h = 2{{k}_{0}}b = 3$ и радиус его основания ${{k}_{0}}a = 1.5,$ расстояние ${{k}_{0}}d = 2.$ Значения диэлектрических и магнитных проницаемостей сред были такие же, что и выше. При решении задачи методом на основе ММОИУ в качестве криволинейных координат $(u,v)$ использовались полярные (u – угловая координата) или цилиндрические координаты ($\rho = v,\,\,z = u$) в случае дифракции на сфере или цилиндре соответственно. Число точек коллокации было выбрано равным ${{N}_{1}} = 24,$ ${{N}_{2}} = 8$ в случае дифракции на сфере и ${{N}_{1}} = 16,$ ${{N}_{2}} = 8$ в случае дифракции на цилиндре. Параметр ${{k}_{0}}\delta = {{10}^{{ - 3}}}.$ Диаграмму рассеяния вычисляли для угла наблюдения $\varphi = 0.$ Как видно из таблицы, результаты расчета диаграммы, полученные методами, основанными на ММДИ и ММОИУ, различаются не более чем на 1 и 1.5% соответственно в случае дифракции на сфере и цилиндре.

Представляет интерес рассмотреть задачу дифракции на неоднородном теле вращения, расположенном в верхнем полупространстве (плоская волна также падает из верхнего полупространства). При этом обе методики решения задачи остаются практически без изменений. Небольшое отличие возникает только при нахождении асимптотики рассеянного поля в верхнем полупространстве. На рис. 3, 4 изображены угловые зависимости модуля диаграммы рассеяния для случая дифракции плоской волны вида (3) (рис. 3а, 3б) или (4) (рис. 4а, 4б). Рассматривалась дифракция на круговом цилиндре, расположенном в верхнем полупространстве. Рис. 3а и 4а соответствуют нормальному падению плоской волны, а рис. 3б и 4б – падению плоской волны под углом ${{\chi }_{1}} = 45^\circ .$ Магнитная проницаемость всюду равна единице, а диэлектрические проницаемости верхнего и нижнего полупространств – ${{\varepsilon }_{1}} = 1,\,\,{{\varepsilon }_{2}} = 4$ соответственно, высота цилиндра $h = {\lambda \mathord{\left/ {\vphantom {\lambda 2}} \right. \kern-0em} 2}$ и радиус его основания $a = {\lambda \mathord{\left/ {\vphantom {\lambda 2}} \right. \kern-0em} 2}.$ Расстояние от центра тела до границы раздела сред $d = {\lambda \mathord{\left/ {\vphantom {\lambda 2}} \right. \kern-0em} 2}.$ Кривые на рисунках построены для полуплоскостей $\varphi = 0$ и $\varphi = \pi .$ Кривые 1 на всех рисунках иллюстрируют угловую зависимость диаграммы рассеяния в верхнем полупространстве при дифракции на неоднородном цилиндре, диэлектрическая проницаемость которого изменяется по закону

(66)
${{\varepsilon }_{i}}(z) = {{\varepsilon }_{c}}\left( {1 + \nu {{{\left( {\frac{z}{b}} \right)}}^{2}}} \right),$
Рис. 3.

Угловая зависимость модуля диаграммы рассеяния цилиндра, расположенного в верхнем полупространстве. Падает плоская волна (3). Нормальное (а) и наклонное (б) падение волны.

Рис. 4.

Угловая зависимость модуля диаграммы рассеяния цилиндра, расположенного в верхнем полупространстве. Падает плоская волна (4). Нормальное (а) и наклонное (б) падение волны.

где εc = 1, ν = 3, b – как и выше, половина высоты. Кривые 2 соответствуют рассеянию на однородном цилиндре с диэлектрической проницаемостью ${{\varepsilon }_{i}} = 2.5.$ Таким образом, диэлектрическая проницаемость однородного цилиндра равна среднему арифметическому значений проницаемости неоднородного тела, взятых при $z = 0$ и $z = b.$ Задачу дифракции на неоднородном цилиндре решали с помощью ММОИУ. Как видно из рисунков, в случае нормального падения волны зависимости диаграммы рассеяния неоднородного и однородного тел близки между собой. В случае наклонного падения волны типа (3) отличие зависимостей друг от друга существенно лишь в окрестности главного лепестка диаграммы рассеяния. Если же структура облучается плоской волной вида (4) и угол ${{\chi }_{1}} = 45^\circ ,$ то зависимости неоднородного и однородного тела совершенно различны во всем диапазоне углов наблюдения.

ЗАКЛЮЧЕНИЕ

Предложены два подхода к решению векторной задачи дифракции на теле вращения, расположенном в однородном диэлектрическом полупространстве. Для проверки работы методики, основанной на ММДИ, построены графики невязки краевого условия на контуре осевого сечения однородной сферы и суперэллипсоида вращения, погруженных в полупространство. Показано, что максимальный уровень невязки для тела, граница которого имеет участки резкого изменения координат (суперэллипсоид), не превосходит $5 \times {{10}^{{ - 4}}}.$ Проведено сравнение методов на примере задачи дифракции на однородной сфере и круговом цилиндре, расположенных в полупространстве. Продемонстрировано, что результаты расчета диаграммы рассеяния, полученные методами, основанными на ММДИ и ММОИУ, различаются не более чем на 1 и 1.5% соответственно в случае дифракции на сфере и цилиндре. Проведены расчеты диаграммы рассеяния для задачи дифракции на круговом цилиндре, расположенном в верхнем полупространстве. Результаты представлены как для однородного цилиндра, так и для неоднородного цилиндра (тех же размеров), диэлектрическая проницаемость которого зависит от вертикальной координаты по квадратичному закону (в цилиндрической системе координат). Диэлектрическая проницаемость однородного тела выбиралась равной среднему значению проницаемости неоднородного цилиндра. Рассмотрены две поляризации падающего поля. Показано существенное отличие диаграмм рассеяния неоднородного и однородного цилиндров для случая наклонного падения плоской волны, у которой вектор электрического поля лежит в плоскости падения.

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

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

  1. Жук Н.П., Яровой А.Г. // ЖТФ. 1992. Т. 62. № 7. С. 1.

  2. Zhuck N.P., Yarovoy A.G. // IEEE Trans. 1994. V. AP-42. № 1. P. 16.

  3. Жук Н.П., Шульга С.Н., Яровой А.Г. // ЖТФ. 1998. Т. 68. № 1. С. 1.

  4. Калинченко Г.А., Кюркчан А.Г., Лерер А.М. и др. // РЭ. 2001. Т. 46. № 9. С. 1087.

  5. Еремин Ю.А., Орлов Н.В. // Оптика и спектроскопия. 1998. Т. 84. № 4. С. 625.

  6. Гришина Н.В., Еремин Ю.А. // Оптика и спектроскопия.1999. Т. 86. № 3. С. 475.

  7. Гришина Н.В., Еремин Ю.А. // Оптика и спектроскопия. 2000. Т. 88. № 2. С. 284.

  8. Кюpкчан А.Г., Маненков С.А., Негорожина Е.С. // РЭ. 2006. Т. 51. № 11. С. 1285.

  9. Anyutine A.P., Kyurkchan A.G., Manenkov S.A., Minaev S.A. // J. Quantitative Spectroscopy and Radiative Transfer. 2006. V. 100. P. 26.

  10. Кюpкчан А.Г., Маненков С.А., Негорожина Е.С. // РЭ. 2008. Т. 53. № 3. С. 276.

  11. Kyurkchan A.G., Manenkov S.A. // J. Quantitative Spectroscopy & Radiative Transfer. 2008. V. 109. P. 1430.

  12. Maнeнкoв C.A. // PЭ. 2008. T. 53. № 7. C. 789.

  13. Maнeнкoв C.A. // PЭ. 2018. T. 63. № 1. C. 3.

  14. Кюркчан А.Г., Смирнова Н.И. Математическое моделирование в теории дифракции с использованием априорной информации об аналитических свойствах решения. М.: Медиа Паблишер, 2014.

  15. Kyurkchan A.G., Manenkov S.A. // J. Quantitative Spectroscopy & Radiative Transfer. 2012. V. 113. P. 2368.

  16. Дмитриев В.И., Захаров Е.В. Интегральные уравнения в краевых задачах электродинамики. М.: Изд-во МГУ, 1987.

  17. Захаров Е.В., Пименов Ю.В. Численный анализ дифракции радиоволн. М.: Радио и связь, 1982.

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