Применение тонкопроволочного интегрального представления электромагнитного поля к решению задачи дифракции электромагнитных волн на проводящих телах

Полный текст

Введение

Термин «дифракция» дословно означает «отклонение». В электродинамике в «широком» смысле слова под дифракцией принято понимать все явления, связанные с распространением электромагнитных полей, создаваемых когерентными источниками, при наличии препятствий различных электрических размеров [1]. Благодаря явлению дифракции электромагнитные волны могут попадать в область геометрической тени, огибать препятствия, стлаться вдоль поверхности, проникать через малые отверстия в экранах и т. д.

Интерес к подобным задачам возник довольно давно. Традиционная теория дифракции создавалась на протяжении нескольких столетий Х. Гюйгенсом, О. Френелем, Г. Гельмгольцем, Г.Р. Кирхгофом и другими авторами. Для понимания волновых процессов и расчета дифракционных полей большое значение имеет принцип Гюйгенса, согласно которому распространение волн обусловлено действием вторичных источников. Френель уточнил принцип Гюйгенса, приняв во внимание интерференцию сферических волн, излучаемых вторичными источниками. Дальнейшее уточнение принципа Гюйгенса – Френеля принадлежит Кирхгофу, который дал его строгую формулировку, основываясь на уравнении Гельмгольца. В случае, когда тело, на которое падает электромагнитная волна, обладает бесконечно большой проводимостью, строгое решение задачи дифракции заключается в ее сведении к векторному интегральному уравнению (ИУ) на поверхности тела. Этот подход часто называют методом поверхностных токов. Здесь можно выделить два этапа. Первый заключается в вычислении распределения поверхностных токов (внутренняя задача электродинамики), а второй – в вычислении поля рассеяния, создаваемого поверхностными токами. Решение внутренней задачи представляет собой достаточно серьезную проблему, в наиболее основательном случае требующую исследований разрешимости ИУ, корректного выбора пространств для решений и т. д. [2].

В современных условиях для решения задачи дифракции используются системы автоматизированного проектирования, такие как CST STUDIO, HFSS и FEKO. Они используют разные методы, среди которых можно отметить метод моментов [3]. С помощью метода моментов исходное векторное ИУ сводится к СЛАУ относительно неизвестных коэффициентов разложения токовой функции в ряд по заранее выбранной системе базисных функций. В наиболее популярном на сегодняшний день варианте поверхность тела подвергается процедуре триангуляции, позволяющей затем использовать в качестве базисных функции Рао – Уилтона – Глиссона [4]. Определение коэффициентов матрицы СЛАУ при этом в общем случае предполагает вычисление интегралов четвертой степени кратности в локальных координатах треугольников, что представляет собой довольно сложную в вычислительном плане задачу. Особым образом производится вычисление диагональных элементов матрицы СЛАУ.

Решение задачи дифракции можно существенно упростить, применяя тонкопроволочное интегральное представление электромагнитного поля (ИП ЭМП) [5]. В этом случае следует использовать ранее осуществленную триангуляцию, заменив грани треугольников на тонкие проводники малого электрического радиуса. При малых электрических размерах треугольников сетка проводников будет эквивалентна сплошной металлической поверхности, вследствие чего можно ожидать хорошей аппроксимации решения исходной задачи дифракции. Внутренняя задача в данном случае с помощью метода коллокаций также сводится к СЛАУ, нахождение элементов матрицы которой заключается в вычислении довольно простых одномерных интегралов. Здесь также следует отметить, что сеточные структуры – это отдельный класс электродинамических структур, а решение для них задач дифракции и излучения имеет очень важное теоретическое и прикладное значение. Ранее в [6; 7] было показано, что с помощью тонкопроволочного ИП ЭМП можно успешно решать задачи дифракции на некоторых классах метаструктур [8; 9].

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

1. Основные расчетные выражения

В наиболее общем виде ЭМП в точке наблюдения $\mathbf{r}$ создаваемое идеально проводящим телом V, ограниченным поверхностью S (рис. 1) и находящемся в однородной изотропной среде, можно описать интегральным представлением следующего вида [1]:

$\begin{array}{l}\mathbf{E}\mathbf{\left(}\mathbf{r}\mathbf{\right)}=\frac{W_m}{ik}{\int }_S\left(k^2+\nabla \nabla \text{}\cdot \text{}\right)G\mathbf{J}\mathbf{(}\mathbf{r}\mathbf{\text{'}}\mathbf{)}d\mathbf{r}\mathbf{\text{'}};\\ \mathbf{H}\mathbf{\left(}\mathbf{r}\mathbf{\right)}=\nabla \text{}\times \text{}{\int }_{S}G\mathbf{J}(r\text{'})d\mathbf{r}\mathbf{\text{'}}.\end{array}$ (1)

 

Рис. 1. Общая постановка задачи

 

Здесь $W_m$ и $k$ – волновое сопротивление и волновое число среды соответственно; $\mathbf{J}\mathbf{(}\mathbf{r}\mathbf{\text{'}}\mathbf{)}$ – поверхностная плотность электрического тока на $S\supset \mathbf{r}\mathbf{\text{'}}$; $\mathbf{r}\mathbf{\text{'}}$ – точка источника; $\mathbf{\nabla }$ – оператор набла, применяющийся к точке наблюдения,

$G=G\left(R\right)=\frac{1}{4\pi }\frac{e^{-ikR}}{R},R=|\mathbf{r}\mathbf{-}\mathbf{r}\mathbf{\text{'}}|$ (2)

– функция Грина для свободного пространства. Решение внутренней задачи подразумевает определение функции$\mathbf{J}$  при заданном распределении сторонних источников электрического поля ${\mathbf{E}}^{\left(\mathrm{ext}\right)}$ на S. Сформулировать внутреннюю задачу в общем случае можно с помощью граничного условия для идеально проводящей поверхности S:

$\forall \mathbf{r}\mathbf{\text{'}}\in S:\tau \cdot ({\mathbf{E}}^{\left(ext\right)}+\mathbf{E})=\mathrm{0,}$ (3)

здесь $\tau$ – вектор касательной к $S$ в точке $\mathbf{r}\mathbf{\text{'}}$.

Отметим, что аналитически решить внутреннюю задачу можно только в ряде наиболее простых случаев, поэтому на практике используют метод моментов [3], в рамках которого неизвестная токовая функция $\mathbf{J}\mathbf{\left(}\mathbf{r}\mathbf{\right)}$ представляется в виде конечного ряда:

$\mathbf{J}\mathbf{\left(}\mathbf{r}\mathbf{\right)}=\sum _{n=1}^N{a}_n{\mathbf{f}}_n\left(\mathbf{r}\right).$ (4)

Здесь ${\mathbf{f}}_n\left(\mathbf{r}\right)$ – известные векторные базисные функции; $a_n$ – неизвестные коэффициенты, подлежащие определению в процессе решения внутренней задачи.

После подстановки (4) в первое выражение (1), скалярного умножения обеих частей на ${\mathbf{f}}_m\left(\mathbf{r}\right)$ переноса операции дивергенции поля с точки наблюдения на точку источника и интегрирования по точкам источника получаем формулировку внутренней задачи в виде СЛАУ относительно неизвестных коэффициентов $a_n$:

$\sum _n^N{z}_{m,n}{a}_n=b_m,m=1\dots N.$ (5)

Здесь

$\begin{array}{l}{z}_{m,n}={\int }_S{\int }_S\left({\mathbf{f}}_m\left(\mathbf{r}\right)\cdot {\mathbf{f}}_n(\mathbf{r}\text{'})k^2-\right.\\ \left.-[\nabla \cdot {\mathbf{f}}_m(\mathbf{r}\left)\right][\nabla \text{'}\cdot {\mathbf{f}}_n(\mathbf{r}\text{'}\left)\right]\right)Gd\mathbf{r}\text{'}d\mathbf{r};\\ b_m=-\frac{W_m}{ik}{\int }_S{\mathbf{f}}_m\left(\mathbf{r}\right)\text{}\cdot \text{}{E}^{\left(ext\right)}\left(\mathbf{r}\right)d\mathbf{r}\end{array}$ (6)

– матричные коэффициенты и коэффициенты правой части СЛАУ соответственно, оператор $\nabla \text{'}$ применяется к точке источника. Как правило, в современных САПР сначала осуществляется триангуляция поверхности S, а затем в качестве базисных используются функции Рао – Уилтона – Глиссона (RWG), носителем которых является пара треугольных элементов, имеющих общую грань. Следует отметить, что определение элементов $z_{m,n}$ в этом случае представляет собой довольно сложную вычислительную процедуру.

Альтернативный подход к решению задачи дифракции на проводящем теле $V$ можно реализовать на основе тонкопроволочного ИП ЭМП. Для этого необходимо использовать полученную ранее триангуляцию исходной поверхности S, при которой образуется $N_B$ граней. Каждую грань можно заменить на тонкий проводник $L_j$ радиуса $ϵ\ll \lambda$, имеющий соответствующую его обозначению длину $(j=1...N_B)$. Под действием стороннего поля на каждом проводнике возникнет распределение полного тока $I_j\left(l\right)$ ЭМП, создаваемое такой структурой можно записать в виде [5]:

$\mathbf{F}\mathbf{\left(}\mathbf{r}\mathbf{\right)}=\sum _{j=1}^{N_B}{\int }_{L_j}{I}_j\left(l^{\text{'}}\right){\mathbf{K}}^{\left(F\right)}(\mathbf{r}\mathbf{,}{\mathbf{r}}_j(l^{\text{'}}\left)\right)d{l}^{\text{'}},F\equiv E,H;$ (7)

здесь $I_j(l\text{'})$ – распределение полного тока по образующей $L_j$,

$\begin{array}{l}{\mathbf{K}}^{\left(E\right)}(\mathbf{r},\mathbf{r}\text{'})=\frac{W_m}{ik}\left[k^2\widehat{\mathbf{l}}\text{'}G(\mathbf{r},\mathbf{r}\text{'})dl-\frac{\partial }{\partial l}\left((\mathbf{r}-\mathbf{r}\text{'})B(\mathbf{r},\mathbf{r}\text{'})\right)\right];\\ {\mathbf{K}}^{\left(H\right)}(\mathbf{r},\mathbf{r}\text{'})=\widehat{\mathbf{l}}\text{'}\times (\mathbf{r}-\mathbf{r}\text{'})B(\mathbf{r},\mathbf{r}\text{'}).\end{array}$ (8)

– ядра ИП ЭМП, $\mathbf{r}\mathbf{\text{'}}\mathbf{=}{\mathbf{r}}_{\mathit{j}}\mathit{\left(}{l}^{\mathit{\text{'}}}\mathit{\right)}$ – векторное уравнение образующей $L_j$; $\widehat{\mathbf{l}}\text{'}={\widehat{\mathbf{l}}}_j\left(l^{\text{'}}\right)=d{\mathbf{r}}_j\left(l^{\text{'}}\right)/d{l}^{\text{'}}$ – единичный вектор касательной, определенный в точке $l\text{'}$ на образующей $L_j$;

$B=\frac{1}{R}\frac{\partial G}{\partial R}=-\frac{ikR+1}{R^2}G$

– производная функции Грина по R;

$R=\sqrt{|\mathbf{r}\mathbf{-}\mathbf{r}\text{'}{|}^2+{ϵ}^2}$

– расстояние, регуляризированное радиусом проводников $\epsilon$.

Представим каждый проводник в виде совокупности $D+1$ узлов:

$L_j^{\left(D\right)}: {\mathbf{r}}_{j,1},{\mathbf{r}}_{j,2}, ..., {\mathbf{r}}_{i,D+1}$

Уравнение сегмента, соединяющего соседние узлы с номерами $m$ и $m+1$ может быть записано в следующем виде:

${\mathbf{r}}_{i,m}^{\left(s\right)}\left(l\right)={\mathbf{r}}_{i,m}^{*}+{\hat{\mathbf{I}}}_{i,m}l, l\in [-{∆}_{i,m}/2]$

Здесь ${\mathbf{r}}_{i,m}^{*}=({\mathbf{r}}_{\mathit{i}\mathit{,}\mathit{m}}\mathbf{+}{\mathbf{r}}_{\mathit{i}\mathit{,}\mathit{m}\mathit{+}\mathit{1}}\mathit{)}\mathit{/}\mathit{2}$ – центр сегмента; ${∆}_{i,m}=\left|{\mathbf{r}}_{i,m+1}-{\mathbf{r}}_{i,m}\right|$ – длина сегмента, ${\hat{\mathbf{I}}}_{i,m}=\left({\mathbf{r}}_{i,m+1}-{\mathbf{r}}_{i,m}\right)/{∆}_{i,m}$ – единичный вектор касательной на сегменте. Распределение тока на каждом сегменте при $∆\ll \lambda$ можно считать равномерным: $I_{i,m}\left(l\right)=I_{i,m}$. После сегментации всех $L_j$ можно ввести сквозной индекс $k$ для сегментов и переписать (7) в виде

$\mathbf{F}\mathbf{\left(}\mathbf{r}\mathbf{\right)}=\sum _{k=1}^{N_s}{I}_k{\int }_{{\mathrm{\Delta }}_{\mathrm{k}}}{\mathbf{K}}^{\left(F\right)}(\mathbf{r}\mathbf{,}{\mathbf{r}}_k^{\left(s\right)}(l^{\text{'}}\left)\right)d{l}^{\text{'}},F\equiv E,H;$ (9)

Применяя граничное условие вида (3) в центре каждого сегмента, получаем СЛАУ для вычисления неизвестных амплитуд токов $I_k$:

$\sum _{k=1}^{N_s}{z}_{p,k}{I}_k=E_k,p=1\dots N_s,$ (10)

в которой:

$\begin{array}{l}{z}_{p,k}={\widehat{\mathbf{l}}}_p\cdot {\int }_{{\Delta }_k}{\mathbf{K}}^{\left(E\right)}({\mathbf{r}}_p^{*},{\mathbf{r}}_k^{\left(s\right)}(l^{\text{'}}\left)\right)d{l}^{\text{'}},\\ E_k=-\frac{W_m}{ik}{\widehat{\mathbf{l}}}_p\cdot {\mathbf{E}}^{\left(ext\right)}\left({\mathbf{r}}_p^{*}\right).\end{array}$ (11)

Корректное и устойчивое решение СЛАУ в рамках метода коллокаций достигается при выполнении условия:

$2\epsilon \le ∆\le 12\epsilon$ (12)

для любого сегмента [10]. Вне всякого сомнения, решение (11) гораздо проще в численном плане, чем решение (6). При этом интересным моментом является сравнение результатов, полученных двумя описанными методами.

2. Исследуемые переизлучающие структуры

В качестве тестовых структур рассмотрим три объекта сферической формы (рис. 2), центр $O$ которых совпадает с центром декартовой системы координат. Объекты образованы тонкими металлическими проводниками радиуса $L_j$ $\epsilon \ll \lambda$. Радиус сферы, описанной вокруг рассматриваемых объектов, обозначим как r. Возбуждение структур будем осуществлять плоской электромагнитной волной (ПЭМВ), распространяющейся вдоль оси  и поляризованной в направлении оси Oz, имеющей единичную амплитуду и нулевую начальную фазу. Совокупность проводников $L_j$ образует поверхность-сетку с ячейками треугольной формы, имеющими незначительные отличия по форме от равностороннего треугольника, длину грани которого мы обозначим как d. Будем считать, что каждый проводник можно представить в виде $D$ сегментов s равной длины $∆\ll \lambda$, а распределение полного тока на каждом из них можно считать равномерным ($\forall s:I\left(l\right)=I$, здесь $l$ – продольная координата на s). Таким образом, от совокупности проводников мы переходим к совокупности сегментов, ЭМП которых описывается выражением (9), а решение внутренней задачи – выражением (11).

 

Рис. 2. Геометрия исследуемых переизлучающих структур

 

В случае когда $d\ll \lambda$ рассматриваемая сетка становится подобна сплошной металлической поверхности, образованной треугольными элементами $T_k$, на которых можно ввести распределения плотности эквивалентных поверхностных электрических токов ${\mathbf{J}}_k\left(\mathbf{r}\right); (\mathbf{r}\in T_k)$. В этом случае для представления поверхностных токов можно использовать RWG-базис, ЭМП структуры будет определяться выражением (1), а решение внутренней задачи – выражением (6).

Здесь также следует отметить, что в работе Дж. Ми [11] представлено строгое решение задачи дифракции ПЭМВ на однородной сфере произвольного диаметра и состава, находящейся в однородной среде.

3. Результаты численного моделирования

Параметры численного моделирования удобно представить в виде таблицы.

 

Таблица. Параметры численного моделирования

Table. Numerical simulation parameters

№ структуры	D	$∆/r$	$\epsilon /r$	$N_s$
1	1	0,075	0,01	7680
2	4	0,075	0,01	1920
3	8	0,073	0,01	960

 

На рис. 3 показано сравнение результатов расчета нормированных диаграмм рассеяния в меридианной $(\phi =0)$ и азимутальной плоскостях для первой структуры при различных соотношениях $r/\lambda$. Графики типа  получены на основе тонкопроволочного ИП ЭМП (ТП-метод), графики типа $b$ – на основе решения с помощью эквивалентных токов, представленных в RWG-базисе (RWG-метод). Как видно из представленных рисунков, полученные результаты практически не имеют визуальных отличий. Небольшие отличия можно заметить лишь при максимальном отношении $r/\lambda$ для значений $\theta$ и $\phi$ соответствующих боковым лепесткам излучения.

 

Рис. 3. Сравнение нормированных ДР в меридианной и азимутальной плоскостях: а – ТП-метод; б – МКЭ; RWG-метод; 1 – $r=0,25\lambda$; 2 – $r=0,5\lambda$; 3 – $r=1\lambda$; 4 – $r=2\lambda$

 

На рис. 4 показано сравнение результатов расчета нормированных ДР в меридинанной и азимутальной плоскостях для структур 1–3, полученных при различных соотношениях $r/\lambda$. Представленные результаты позволяют оценить влияние размеров ячейки на характеристики излучения структуры. Из рисунка видно, что наибольшие отличия наблюдаются в меридианной плоскости в направлениях, отличных от направления главного лепестка. При этом при увеличении размеров ячейки увеличивается число лепестков бокового излучения, а с ростом частоты их распределение становится асимметричным, что связано с погрешностями в геометрии рассматриваемых структур.

 

Рис. 4. Нормированные меридианные (слева) и азимутальные (справа) ДН для структур 1–3 (номер кривой соответствует номеру структуры): а – $r=0,25\lambda$; б – $r=0,5\lambda$; в – $r=1,0\lambda$

 

Заключение

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

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

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

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

Об авторах

Дмитрий Петрович Табаков

Поволжский государственный университет телекоммуникаций и информатики

Автор, ответственный за переписку.
Email: illuminator84@yandex.ru

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

Область научных интересов: электродинамика, устройства СВЧ и антенны, оптика, численные методы математического моделирования.

Россия, 443010, г. Самара, ул. Л. Толстого, 23

Сергей Владимирович Морозов

Поволжский государственный университет телекоммуникаций и информатики

Email: grimmxxx@gmail.com

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

Область научных интересов: электродинамика, устройства СВЧ и антенны, оптика, численные методы математического моделирования.

Россия, 443010, г. Самара, ул. Л. Толстого, 23

Дмитрий Сергеевич Клюев

Поволжский государственный университет телекоммуникаций и информатики

Email: klyuevd@yandex.ru
ORCID iD: 0000-0002-9125-7076

доктор физико-математических наук, доцент, заведующий кафедрой радиоэлектронных систем Поволжского государственного университета телекоммуникаций и информатики, г. Самара, Россия. Автор более 200 научных работ.

Область научных интересов: электродинамика, устройства СВЧ, антенны, метаматериалы.

Россия, 443010, г. Самара, ул. Л. Толстого, 23

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

Дополнительные файлы