Решение внутренней задачи для конечной регулярной двумерной решетки спиральных элементов, возбуждаемой плоской электромагнитной волной

Полный текст

Введение

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

Актуальность обозначенной задачи, обозначенной в самом начале, связана с ее очень большой общностью. Так, под многоэлементными излучающими структурами можно понимать любые структуры, имеющие более одного элемента – носителя токовых функций. Под это определение попадает подавляющее большинство используемых на практике антенн, и наиболее ярким примером здесь являются фазированные антенные решетки. Под определение многоэлементных переизлучающих структур попадают киральные структуры [1], фотонные кристаллы [2], дифракционные решетки, некоторые типы замедляющих структур и линз диапазона СВЧ, а также любые материальные среды, рассматриваемые на микроскопическом уровне. Все перечисленные переизлучающие структуры, за исключением последних можно отнести к метаструктурам (метаматериалам). Метаматериалы синтезируются внедрением в исходный природный материал частиц с самыми различными формами, которые модифицируют диэлектрическую и магнитную восприимчивость исходного материала. Свойства метаматериалов обусловлены периодической микроструктурой, создаваемой человеком. Основополагающей работой по метаматериалам можно считать работу [3] В.Г. Веселаго. Метаматериалы могут обладать свойствами, принципиально отсутствующими у естественных материалов, встречающихся в природе, например – отрицательный коэффициент преломления [4] (левосторонние материалы). Подобные материалы были изобретены и продемонстрированы Д.Р. Смитом и сотрудниками калифорнийского университета в Сан-Диего [5] спустя более чем 30 лет после публикации статьи В.Г. Веселаго. В основу были положены работа Дж.В. Пендри [6–8]. В настоящее время метаматериалы находят применение в СВЧ- и антенно-фидерных устройствах, позволяя улучшить их параметры. В обзоре [9], представлены результаты, касающиеся применения метаматериалов при построении СВЧ-устройств и антенн. Примером метаматериалов другого рода являются фотонные кристаллы [2], на основе которых можно создавать оптические фильтры, волноводы и многие другие устройства.

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

На протяжении развития электродинамики, оптики и техники СВЧ-устройств появилось достаточно большое количество соответствующих методов, имеющих свои достоинства и недостатки. Долгое время в инженерной практике для расчета взаимодействия в многоэлементных излучающих структурах использовался метод наведенных ЭДС. Основные принципы этого метода изложены в работе [10]. Метод наведенных ЭДС позволяет находить наведенные собственные сопротивления элементов антенны, а также амплитуды и фазы токов в пассивных элементах. Из недостатков метода отметим, что он накладывает определенные ограничения на длину и расстояние между элементами в силу использования приближенных распределений тока. Сейчас для решения подобных задач применяют системы автоматизированного проектирования (САПР), использующие в своей основе метод моментов [11], метод конечных элементов и метод конечных разностей [12]. К недостаткам такого подхода можно отнести высокие требования к ЭВМ, высокую стоимость программного обеспечения, отсутствие математической модели анализируемой структуры в явном виде.

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

Таким образом, разработка универсальных методов расчета межэлементного взаимодействия является актуальной задачей. В [14] рассмотрен итерационный подход к решению интегральных уравнений теории проволочных антенн на основе многошагового метода минимальных невязок. Здесь метод применялся непосредственно к общей матрице системы линейных алгебраических уравнений (СЛАУ). В [15] в качестве основы для расчета взаимодействия было предложено использование модификации метода Гаусса – Зейделя [16] для случая блочной матрицы СЛАУ. Метод применен для решения задачи дифракции на метаструктуре - тонком слое конечных размеров, состоящем из двойных разомкнутых колец. Сделан вывод о том, что предложенный метод может быть эффективен для расчета метаструктур, имеющих конечные размеры. В [17]: рассмотрена задача дифракции плоской электромагнитной волны на киральном слое, представляющим собой прямоугольную эквидистантную решетку 10 × 10 элементов, состоящую из пятисегментных S-элементов. Показано, что несмотря на довольно малое количество хаотически ориентированных S-элементов, диаграммы рассеяния имеют вполне детерминированный характер. В [18] материалы статьи [15] были дополнены алгоритмами расчета элементов блочной матрицы СЛАУ для случая построения структуры из однотипных элементов, обладающих различными видами симметрий. Показано, что данные алгоритмы позволяют существенно сокращать время расчета. В [19] рассмотрен вопрос использования метода [15] для расчета распределений токов на директорной антенне [20]. Показано, что для директорной антенны в исследуемом диапазоне частот метод обладает хорошей сходимостью. В [21] исследована возможность использования итерационного подхода к решению внутренней задачи для симметричного вибратора с рефлектором из параллельных прямолинейных проводников. Вместе с этим представлен ключевой алгоритм расчета элементов блочной матрицы, имеющий линейную оценку сложности в зависимости от числа элементов рефлектора, даны рекомендации по выбору систем проекционных функций и учету имеющихся у структуры симметрий. Проведено исследование распределений токов, входного сопротивления и характеристик излучения рассматриваемой структуры.

В настоящей статье рассмотрено решение внутренней задачи для конечной регулярной двумерной решетки, состоящей из конических спиральных элементов. Возбуждение решетки осуществлялось плоской линейно-поляризованной электромагнитной волной. Решение внутренней задачи осуществлялось на основе метода, предложенного в [15]. Предложен алгоритм вычисления компактной блочной матрицы на основе алгоритмов, описанных в [18]. С учетом влияния резонансных свойств элементов, образующих решетку, на сходимость итерационного процесса, а также с целью обоснования выбора системы проекционных функций в рамках метода моментов, проведен спектральный анализ интегрального оператора внутренней задачи для одиночного спирального элемента, суть которого описана в [22]. Приведены результаты спектрального анализа токовых функций для одиночного спирального элемента и для элементов, находящихся в решетке. Исследована сходимость решения внутренней задачи для решетки элементов. Полученным результатам дана соответствующая интерпретация. Определены направления дальнейшего развития темы.

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

Основные расчетные выражения были достаточно подробно изложены в [21]. Здесь мы изложим соответствующий материал с уточнением некоторых моментов. Рассмотрим излучающую либо переизлучающую структуру, образованную совокупностью элементов $v_n,$ $(n\in \mathcal{N}\text{:}1\dots N, N$  – число элементов), полученных на основе совокупности уникальных элементов ${\underset{¯}{v}}_u$ путем параллельного переноса и (или) вращения $(u\in \mathcal{U}\text{:}1\dots U, U$  – число уникальных элементов). Внутренняя задача электродинамики для рассматриваемой структуры сводится к операторной системе вида [19]:

$\begin{array}{l}\sum _{n^{\text{'}}}{\Psi }_{n\mathrm{,}{n}^{\text{'}}}{\eta }_{n^{\text{'}}}{\zeta }_n\mathrm{,}{n}^{\text{'}}\in \mathcal{N}\mathrm{.}\end{array}$   (1)

В данной системе ${\Psi }_{n\mathrm{,}{n}^{\text{'}}},$ $n\ne n^{\text{'}}$ – операторы взаимодействия элементов $v_n$ и $v_{n^{\text{'}}};$ ${\Psi }_{n\mathrm{,}n}={\Psi }_n$ – собственный оператор элемента $v_n;$ ${\zeta }_n$ – функции, порождаемые сторонними полями, ${\eta }_{n^{\text{'}}}$ – токовые функции, подлежащие определению. К системе (1) применяется схема метода моментов [11]: с помощью систем базисных ${\left\{\beta \right\}}_{n^{\text{'}}}\equiv {{\beta }_{q^{\text{'}}}}_{n^{\text{'}}}$ и тестовых ${\tau }_n\equiv {{\tau }_q}_n$ функций $(q=q_n\in \mathcal{Q}\text{: }1\dots Q\mathrm{,} q^{\text{'}}=q_{n^{\text{'}}}\in \mathcal{Q}\text{: }1\dots Q^{\text{'}}),$ определяется оператор скалярного произведения функций. Искомые функции ${\eta }_{n^{\text{'}}}$ аппроксимируются рядами:

${\eta }_{n^{\text{'}}}=\sum _{q^{\text{'}}}{I}_{q^{\text{'}}}^{\left(n^{\text{'}}\right)}{\beta }_{n^{\text{'}},q^{\text{'}}}.$   (2)

Относительно неизвестных коэффициентов $I_{q^{\text{'}}}^{n^{\text{'}}}$ формируется система линейных алгебраических уравнений (СЛАУ) вида:

$\overline{Z}I=E.$   (3)

Здесь $\overline{Z}$ – блочная матрица обобщенных импедансов, $I,$ $E$ – соответственно блочные векторы неизвестных и правой части,

$\overline{Z}\supset {\overline{z}}_{n,n^{\text{'}}}\supset Z_{q,q^{\text{'}}}^{(n,n^{\text{'}})}=<{\tau }_{n,q},{\Psi }_{n,n^{\text{'}}}\left({\beta }_{n^{\text{'}},q^{\text{'}}}\right)>,$

$I\supset i_{n^{\text{'}}}\supset I_{q^{\text{'}}}^{\left(n^{\text{'}}\right)},E\supset e_n\supset E_q^{\left(n\right)}=<{\tau }_{n,q},{\zeta }_n>.$

Элементы блочной матрицы с одинаковыми индексами ${\overline{z}}_{n,n}={\overline{z}}_n$ будем называть матрицами собственных импедансов, элементы с разными индексами – матрицами взаимных импедансов. Последние определяют взаимодействия между элементами метаструктуры. Скобки «$<>$» используются для обозначения скалярного произведения. Введем в рассмотрение массив уникальных матриц $\left\{{\underset{¯}{\overline{z}}}_k\right\}\equiv \left\{\underset{¯}{\overline{z}}\right\} (k\in \mathcal{K}\text{: }1\dots U\dots K).$   Будем считать, что первые $U$ матриц массива $\left\{\underset{¯}{\overline{z}}\right\}$ являются собственными матрицами импедансов уникальных элементов ${\underset{¯}{v}}_u,$ а оставшиеся $K-U$ элементов – матрицами взаимных импедансов. Массив $\left\{p_{n,n^{\text{'}}}^{\left(\overline{z}\right)}\right\}$ с элементами $p_{n,n^{\text{'}}}^{\left(\overline{z}\right)}=k$будет устанавливать соответствие между парой индексов $n\mathrm{,}{n}^{\text{'}}\in \mathcal{N}$ и индексом $k\in \mathcal{K}.$ Таким образом, имеем:

${\overline{z}}_{n,n^{\text{'}}}={\underset{¯}{\overline{z}}}_k,k=p_{n,n^{\text{'}}}^{\left(\overline{z}\right)}\in \left\{p_{n,n^{\text{'}}}^{\left(\overline{z}\right)}\right\}\equiv \left\{p^{\left(\overline{z}\right)}\right\}.$

В общем случае, когда структура состоит из хаотично расположенных, отличных друг от друга элементов $v_n,$ введение представленных массивов $\left\{\underset{\_}{\overline{z}}\right\}$ и $\left\{p^{\left(\overline{z}\right)}\right\}$ не имеет смысла, и лишь несущественно увеличивает затраты машинной памяти и машинного времени. Но в случае, когда в структуре имеются однотипные элементы, а также присутствует регулярность в их взаимном расположении, имеет место неравенство $K<N^2,$ и введение указанных массивов становится оправданным, при этом они, по сути, несут в себе базовую часть априорной информации о структуре. Алгоритм расчета массива уникальных матриц будет рассмотрен далее.

Заключительным этапом решения внутренней задачи является решение блочной СЛАУ (3). В данном случае наиболее рациональным представляется использование итерационных методов, основанных на расщеплении матрицы [16], однако в случае отсутствия диагонального преобладания в матрице $\overline{Z}$ они не обеспечивают сходимости итерационного процесса.

Достичь сходимости классических итерационных методов можно рассматривая именно блочную СЛАУ (3) в предположении, что между элементами структуры отсутствует достаточно сильная связь. Также на сходимость итерационного процесса существенным образом влияет выбор систем проекционных функций (СПФ). Оптимальным вариантом в данном случае представляется использование систем собственных функций (СФ) входящих в структуру элементов или СПФ, близких к СФ [23]. Вычисление СФ представляет собой вспомогательную задачу, сложность решения которой определяется числом уникальных элементов ${\underset{¯}{v}}_u$ и их свойствами. Полная проблема собственных значений для матриц собственных импедансов уникальных элементов ${\underset{¯}{\overline{z}}}_u,$ полученных в рамках метода моментов с помощью базисных функций ${\beta }_{u\mathrm{,}q}$ и тестовых функций ${\tau }_{u\mathrm{,}q},$ формулируется в виде:

${\underset{¯}{\overline{z}}}_u{\overline{J}}^{\left(u\right)}=\mathcal{D}\left(x^{\left(u\right)}\right){\overline{J}}^{\left(u\right)}.$

Здесь ${\overline{J}}^{\left(u\right)}$ – матрица, столбцы которой содержат собственные векторы (СВ) матрицы ${\underset{¯}{\overline{z}}}_u,$ $\mathcal{D}\left(x^{\left(u\right)}\right)$ – диагональная матрица, формируемая вектором $x^{\left(u\right)},$ элементы ${\xi }_q$ которого являются соответствующими собственными значениями (СЗ) матрицы ${\underset{¯}{\overline{z}}}_u; \mathcal{D}$ – оператор, формирующий диагональную матрицу. Операндом $\mathcal{D}$ является вектор элементов главной диагонали, либо матрица, элементы главной диагонали которой используются оператором $\mathcal{D}$ для создания диагональной матрицы. Далее будем считать, что для рассматриваемых одиночных элементов комплексные матрицы ${\underset{¯}{\overline{z}}}_u$ симметричны. СЗ ${\xi }_q$ аппроксимируют СЗ собственных интегральных операторов ${\underset{¯}{\Psi }}_u,$ СФ указанных интегральных операторов аппроксимируются с помощью СВ следующим образом:

${\stackrel{⌣}{\beta }}_{u,q}\approx \sum _{q^{\text{'}}=1}^{Q^{\text{'}}}{J}_{q^{\text{'}},q}^{\left(u\right)}{\beta }_{u,q^{\text{'}}},J_{n,q}^{\left(u\right)}\in {\overline{J}}^{\left(u\right)}.$

Далее эта задача будет рассмотрена для конического спирального элемента, образующего исследуемую структуру.

Применяя процедуру Гаусса – Зейделя [16] к блочной СЛАУ, можно записать следующую форму ее решения:

$i_n^{(l+1)}={\overline{p}}_n{e}_n-\sum _{l>i}{\overline{w}}_{n,n^{\text{'}}}{i}_{n^{\text{'}}}^{\left(l\right)}-\sum _{l<i}{\overline{w}}_{n,n^{\text{'}}}{i}_{n^{\text{'}}}^{(l+1)},$   (4)

здесь и далее $l$ – номер шага итерационного процесса,

${\overline{p}}_n={\overline{z}}_{n,n}^{-1},{\overline{w}}_{n,n^{\text{'}}}={\overline{p}}_n{\overline{z}}_{n,n^{\text{'}}}$

– соответственно обращенные собственные матрицы и весовые матрицы, ${\overline{p}}_n$ выступают в качестве матриц предобуславливателя. В процедуре (4) можно также использовать более простой вариант

${\overline{p}}_n=\mathcal{D}{\left({\overline{z}}_{n,n}\right)}^{-1}.$

Если положить во второй сумме (4) $l+1=l,$ то формула будет соответствовать методу простой итерации. Критерий оценки сходимости строится в соответствии с неравенством:

${\delta }_l=\underset{n}{\max }\left(\right|i_n^{(l+1)}-i_n^{\left(l\right)}|/|i_n^{(l+1)}\left|\right)\le {\delta }_{*},$   (5)

где ${\delta }_{*}$ – сколь угодно малое наперед заданное число. Под $\left|v\right|$ здесь и далее будем понимать евклидову норму для вектора $v.$

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

Многоэлементная тонкопроволочная структура $L$ представляет собой совокупность $N$ тонких проводников $L_1\mathrm{,}$ $L_2\mathrm{,}$ $\dots \mathrm{,}$ $L_N$ произвольной формы, расположенных в свободном пространстве с волновым сопротивлением $W_m.$ Для простоты предположим, что радиус всех проводников одинаков и равен $\epsilon .$ Каждый проводник можно описать векторным уравнением, зависящим от натурального параметра l:

$r_n\left(l\right)=\widehat{x}{X}_n\left(l\right)+\widehat{y}{Y}_n\left(l\right)+\widehat{z}{Z}_n\left(l\right),l\in [L_{n,\min },L_{n,\max }],$

где $X_{n}l$ ,$Y_{n}l$, $Z_{n}l$ – гладкие функции. Под $L_n=L_{n\mathrm{,}\max }-L_{n\mathrm{,}\min }$ также будем понимать длину n-го проводника. Интегральное представление электромагнитного поля (ИП ЭМП) такой структуры можно записать в виде [24]:

$F\left(r\right)=\sum _{n^{\text{'}}=1}^N{\int }_{L_{n^{\text{'}}}}{I}_{n^{\text{'}}}\left(l^{\text{'}}\right)K^{\left(F\right)}(r,r_{n^{\text{'}}}(l^{\text{'}}\left)\right)d{l}^{\text{'}},F\equiv E,H;$   (6)

здесь $I_{n^{\text{'}}}{l}^{\text{'}}$ – распределение полного тока по образующей $L_{n^{\text{'}}},$

$K^{\left(E\right)}=\frac{W_m}{ik}\left[k^2\widehat{l}\text{'}Gdl-\frac{\partial }{\partial l}\left((r-r\text{'})B\right)\right],$

$K^{\left(H\right)}=\widehat{l}\text{'}\times (r-r\text{'})B$

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

$G=\frac{\exp \left(-ikR\right)}{4\pi R},B=\frac{1}{R}\frac{\partial G}{\partial R}=-\frac{ikR+1}{R^2}G,$

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

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

Для $I_{n^{\text{'}}}{l}^{\text{'}}$ целесообразно представление в виде рядов типа (2). Исходное ИП ЭМП (6) при этом приобретет вид:

$F\left(r\right)=\sum _{n^{\text{'}}}\sum _{q^{\text{'}}}{I}_{q^{\text{'}}}^{\left(n^{\text{'}}\right)}{\int }_{L_{n^{\text{'}}}}{\beta }_{n^{\text{'}},q^{\text{'}}}\left(l^{\text{'}}\right)K^{\left(F\right)}(r,r_{n^{\text{'}}}(l^{\text{'}}\left)\right)d{l}^{\text{'}},F\equiv E,H;$

На каждой образующей справедливо граничное условие для идеального проводника:

$\left(E^{\left(\right(in\left)\right)}\right(r_n\left(l\right))+E(r_n\left(l\right)\left)\right)\cdot {\widehat{l}}_n\left(l\right)=0.$  (7)

Умножая поочередно (7) на тестовые функции ${\tau }_{n\mathrm{,}q}l$ и интегрируя по l, получаем СЛАУ для вычисления $I_{q^{\text{'}}}^{n^{\text{'}}},$ по форме совпадающую с (3), в которой:

$\begin{array}{l}{Z}_{q,q^{\text{'}}}^{(n,n^{\text{'}})}=\underset{L_n}{\int }\text{}\underset{L_{n^{\text{'}}}}{\int }{\tau }_{n,q}\left(l\right){\beta }_{n^{\text{'}},q^{\text{'}}}\left(l^{\text{'}}\right){\mu }^{(n,n^{\text{'}})}(l,l^{\text{'}})d{l}^{\text{'}}dl,\end{array}$  (8)

$E_q^{\left(n\right)}=\underset{L_n}{\int }{\tau }_{n,q}\left(l\right){\nu }^{\left(n\right)}\left(l\right)dl;$

${\mu }^{(n,n^{\text{'}})}(l,l^{\text{'}})={\widehat{l}}_n\left(l\right)\cdot K^{\left(E\right)}\left(r_n\right(l),r_{n^{\text{'}}}(l^{\text{'}}\left)\right),$

${\nu }^{\left(n\right)}\left(l\right)={\widehat{l}}_n\left(l\right)\cdot E^{\left(\right(in\left)\right)}\left(r_n\right(l\left)\right).$

Конкретизируем способ вычисления интегралов в (8) с помощью процедуры сегментации проводников, предполагающей представление n-го проводника в виде совокупности $M+1$ узлов $L_n^{\left(M\right)}:$ ${\mathrm{r}}_{\mathrm{n},1}, {\mathrm{r}}_{\mathrm{n},2}, \dots \mathrm{,} {\mathrm{r}}_{\mathrm{n},\mathrm{M}+1}.$ Между узлами с номерами $m$ и $m+1$ находится m-й сегмент $r_{n,m}\left(l\right),$ уравнение которого может быть записано в следующем виде:

$r_{n,m}\left(l\right)={\dot{r}}_{n,m}+{\widehat{l}}_{n,m}l,l\in [-{\Delta }_{n,m}/\mathrm{2,}{\Delta }_{n,m}/2].$

Здесь ${\dot{r}}_{n,m}=(r_{n,m}+r_{n,m+1})/2$ – центр сегмента, ${\Delta }_{n,m}=|r_{n,m+1}-r_{n,m}|$ – длина сегмента, ${\widehat{l}}_{n,m}=(r_{n,m+1}-r_{n,m})/{\Delta }_{n,m}$ – единичный вектор касательной на сегменте. Здесь и далее индексы $m\in \mathcal{M}: 1\dots M$ и $m^{\text{'}}\in \mathcal{M}\text{'}: 1\dots M^{\text{'}}$ будем ассоциировать с сегментами на n-ом и $n^{\text{'}}$-ом элементе структуры соответственно.

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

${\beta }_{n,q}^{\left(M\right)}\left(l\right)=\sum _m{\beta }_{n,q}\left({\dot{l}}_{n,m}\right)\sigma (l^{\text{'}},{\dot{l}}_{n,m},{\Delta }_{n,m}),$   (9)

где ${\dot{l}}_{n\mathrm{,}m}$ – значение натурального параметра на сегментированной образующей $L_n^{\left(M\right)},$ соответствующее центру сегмента с индексом m, $\sigma (l,\dot{l},\Delta )$ – функция, описывающая прямоугольный единичный импульс, имеющий центр в точке $\dot{l}$ и ширину $\Delta .$ В качестве тестовых будем использовать взвешенные суммы дельта-функций Дирака:

${\tau }_{n,q}^{\left(M\right)}\left(l\right)=\sum _m{\tau }_{n,q}\left({\dot{l}}_{n,m}\right)\delta (l-{\dot{l}}_{n,m}).$   (10)

Такой подход можно рассматривать как обобщенный метод коллокаций [25]. Из представленных выражений видно, что в случае (9) роль весовых коэффициентов играют значения функций ${\beta }_{n,q}\left(l\right),$ вычисленные в точках коллокации ${\dot{l}}_{n\mathrm{,}m},$ а в случае (10) аналогичная роль принадлежит функциям ${\tau }_{n,q}\left(l\right).$

Применяя приведенные выражения в (8) с учетом свойств дельта-функции, получаем формулы для расчета матричных коэффициентов и коэффициентов правой части СЛАУ с помощью конечных сумм:

$Z_{q,q^{\text{'}}}^{(n,n^{\text{'}})}\approx \sum _m\sum _{m^{\text{'}}}{\tau }_{n,q}\left({\dot{l}}_{n,m}\right){\beta }_{n^{\text{'}},q^{\text{'}}}\left({\dot{l}}_{n^{\text{'}},m^{\text{'}}}\right){\mu }_{m,m^{\text{'}}}^{(n,n^{\text{'}})},$   (11)

$E_q^{\left(n\right)}\approx \sum _m{\tau }_{n,q}\left({\dot{l}}_{n,m}\right){\nu }_m^{\left(n\right)},$

здесь:

${\mu }_{m,m^{\text{'}}}^{(n,n^{\text{'}})}={\int }_{{\Delta }_{n^{\text{'}},m^{\text{'}}}}\text{}{\mu }^{(n,n^{\text{'}})}({\dot{l}}_{n,m},l^{\text{'}})d{l}^{\text{'}},{\nu }_m^{\left(n\right)}={\nu }^{\left(n\right)}\left({\dot{l}}_{n,m}\right).$

Традиционному методу коллокаций соответствует выбор:

${\beta }_{n,q}\left(l\right)={\delta }_{l,{\dot{l}}_{n,q}},{\tau }_{n,q}\left(l\right)={\delta }_{l,{\dot{l}}_{n,q}},$

где ${\delta }_{x\mathrm{,}y}$ – дельта Кронекера. Корректное решение СЛАУ в рамках метода коллокаций достигается при выполнении условия [26] для любого сегмента:

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

3. Исследуемая структура

Геометрия исследуемой метаструктуры (общий вид) приведена на рис. 1, а, геометрия образующих ее элементов – на рис. 1, б. Структура представляет из себя решетку с прямоугольной сеткой, лежащую в плоскости $xOy.$ Общий размер метаструктуры можно определить как $L_x\times L_y.$ В узлы решетки помещены переизлучающие элементы, имеющие одинаковую форму и одинаковую пространственную ориентацию, что, с учетом фиксированного межэлементного расстояния, позволяет считать решетку регулярной. Размер переизлучающих элементов не превышает значения D. Возбуждение переотраженного излучения осуществляется плоской электромагнитной волной (ПЭМВ), имеющей линейную поляризацию. Угол падения ПЭМВ на структуру может иметь произвольное значение. Стороннее электрическое поле $E^{\left(in\right)},$ порождаемое волной, возбуждает токи, распределяющиеся по переотражающим элементам. Эти токи становятся источниками ЭМП, создаваемого метаструктурой.

 

Рис. 1. Геометрия решетки (а) и образующих ее элементов (б)

Fig. 1. Geometry of the lattice (a) and its constituent elements (b)

 

Проведем детальное рассмотрение решетки метаструктуры. Для узлов решетки и соответствующих элементов решетки удобно использовать двойную индексацию. Индексы вдоль соответствующих осей будем использовать обозначения $i_x,$ $i_y.$ При заполнении итоговой матрицы СЛАУ осуществляется переход от двойной к сквозной индексации $(n=1\dots N):$

$n=i_x+(i_y-1)N_x,i_x=\mathrm{1,}\dots ,N_x,i_y=\mathrm{1,}\dots ,N_y.$   (13)

Рассмотрим параметры, задающие геометрию решетки. В дальнейшем $h_x$ и $h_y$ будут означать расстояние между соседними элементами. С использованием этих обозначений можно вычислить координаты узлов решетки:

$r_{i_x,i_y}=\widehat{x}\left(h_x\right(i_x-1)-L_x/2)+\widehat{y}\left(h_y\right(i_y-1)-L_y/2).$

Обозначения $N_x$ и $N_y$ будем использовать для числа элементов вдоль соответствующих осей. В этом случае:

$L_x=(N_x-1)h_x,L_y=(N_y-1)h_y.$

Элементы $L_{i_x\mathrm{,}{i}_y}$ метаструктуры – конические спирали, образованные идеальным проводником, имеющим круглое сечение. Ось каждой спирали спирали параллельна оси $Oz.$ Сам элемент показан на рис. 1, б. Минимальный радиус спирали – $r_{\min },$ максимальный радиус спирали – $r_{\max }.$ Высота спирали равна $H$ число витков обозначим как $N_l.$ Осью спиральных элементов будем считать линию, совпадающую с осью $Oz.$ Диаметр проводников $2\epsilon$ существенно меньше $\lambda$ и существенно меньше размера спиралей и межвиткового расстояния. Это позволяет использовать тонкопроволочное приближение при построении модели метаструктуры (выражения, приведенные в предыдущем разделе).

Общее уравнение конической архимедовой спирали имеет вид:

$L:r\left(t\right)=(r_{\max }-h_{r}t)\cos \left(st\right)\widehat{x}+$   (14)

$\begin{array}{l}+(r_{\max }-h_{r}t)\sin \left(st\right)\widehat{y}+h_{z}t\widehat{z},\\ t\in \left[\mathrm{0,2}\pi N_l\right].\end{array}$

Здесь $t$ является ненатуральным параметром, $r_{\min },$ $r_{\max }$ – соответственно минимальный и максимальный радиус спирали, $H$ – ее высота, $N_l$ – число витков. С помощью $h_r$ обозначен радиальный коэффициент намотки, с помощью $h_z$ – коэффициент намотки по высоте. Параметр $s=\pm 1$ позволяет изменять направление намотки спирали, что дает возможность использовать приведенное параметрическое уравнение для анализа соответствующих киральных структур (левых и правых). Переменные $h_r$ и $h_z$ можно вычислить при заданных $r_{\min },$ $r_{\max },$ $H$ и $N_l$ с помощью следующих формул:

$h_r=\frac{r_{\max }-r_{\min }}{2\pi N_l},h_z=\frac{H}{2\pi N_l}.$   (15)

Для определения натурального параметра на конической спирали следует использовать выражение:

$l\left(t\right)=\underset{0}{\overset{t}{\int }}\left|l\right(t^{\text{'}}\left)\right|d{t}^{\text{'}},$   (16)

в котором $\widehat{l}\left(t\right)=dr\left(t\right)/dt$ – касательный орт, определенный в точке $r\left(t\right)$ образующей. Таким образом, уравнение образующей конической архимедовой спирали, имеющего запись в натуральном параметре, можно получить с помощью (14), где следует использовать подстановку $t\to t\left(l\right).$ Для функции $t\left(l\right)$ на конической архимедовой спирали не существует явного выражения, поэтому она может быть определена только численно с использованием выражения (16) на основе метода обратной интерполяции [16]. Любой элемент $L_{i_x\mathrm{,}{i}_y}$ решетки может быть получен согласно выражению:

$L_{i_x,i_y}:r_{i_x,i_y}\left(l\right)={\overline{R}}^{\left(z\right)}\left({\varphi }_{i_x,i_y}\right)r\left(l\right)+r_{i_x,i_y}.$   (17)

Здесь ${\overline{R}}^{\left(z\right)}\left(\varphi \right)$ – матрица поворота вокруг оси $Oz$ на угол $\varphi .$ Далее можно перейти к сквозной нумерации (13) и использовать для расчетов выражения, приведенные в предыдущем разделе.

Поле ПЭМВ $W$ (рис. 2) имеет следующий вид:

$E^{\left(in\right)}\left(r\right)=\widehat{p}{E}_0\exp (-ikr+{\psi }_0),$   (18)

$H^{\left(in\right)}\left(r\right)=[\widehat{k}\text{\hspace{0.17em}}\times \text{\hspace{0.17em}}\widehat{p}]\frac{E_0}{W_m}\exp (-ikr+{\psi }_0).$

Здесь $k=\widehat{k}k$, $\widehat{k}$ является волновым ортом, $k$ используется для обозначения волнового числа, $E_0$ – амплитуда вектора электрического поля, $W_m$ – волновое (характеристическое) сопротивление среды, ${\psi }_0$ – начальная фаза волны, $\widehat{p}$ – вектор поляризации, определяющий направление колебаний вектора $E$ волны.

 

Рис. 2. Связь плоской электромагнитной волны $\mathit{W}$ с первичной волной ${\mathit{W}}^{\mathbf{\text{'}}}$ и глобальной декартовой системой координат

Fig. 2. Relationship between a plane electromagnetic wave $\mathit{W}$ and the primary wave ${\mathit{W}}^{\mathbf{\text{'}}}$ and the global Cartesian coordinate system

 

Расчеты ПЭМВ $W$ предполагают ее определение пятью параметрами: амплитудой $E_0,$ фазой ${\psi }_0,$ углом поляризации ${\varpi }_0,$ и углами ${\theta }_0$ и ${\phi }_0.$ Эти углы определяют направление распространения волны. Дадим определение базовой волны $W^{\text{'}}$ (рис. 2). Ее орты $\widehat{k}\text{'}, \widehat{p}\text{'}$ определяются ортами глобальной прямоугольной системы координат. Для определенности будем полагать, что $\widehat{p}\text{'}=\widehat{z}, \widehat{k}\text{'}=\widehat{x}.$  Таким образом, волна $W^{\text{'}}$ обладает вертикальной поляризацией и распространяется вдоль оси $Ox.$ Связь ортов волн $W$ и $W^{\text{'}}$ имеет следующий вид:

$v={\overline{R}}^{\left(z\right)}\left({\phi }_0\right){\overline{R}}^{\left(y\right)}\left({\theta }_0\right){\overline{R}}^{\left(x\right)}\left({\varpi }_0\right)v\text{'},v\equiv \widehat{p},\widehat{k}.$

В данном выражении ${\overline{R}}^{\left(f\right)}\left(\varphi \right)$ – матрицы поворота вокруг соответствующих осей, $\varphi$ – угол поворота; $f\equiv x,y,z.$

4. Алгоритм расчета массива уникальных матриц

В [18] рассмотрены непосредственный и ключевой алгоритмы расчета элементов $\left\{p^{\left(\overline{z}\right)}\right\}$ и $\left\{\underset{¯}{\overline{z}}\right\}.$ Указанные алгоритмы, в свою очередь, могут быть строгими и приближенными. Достоинством непосредственного алгоритма является универсальность, а основным недостатком – повышенная вычислительная сложность, связанная с необходимостью сравнения блоков матрицы $\overline{Z}.$ Однако в [18] не было указано, что подобный подход может быть эффективен при расчете на множестве частот и (или) при большом количестве вариантов возбуждения структуры, так как процедура расчета применяется однократно. Для расчета массива уникальных матриц рассматриваемой структуры предлагается использование приближенного комбинированного алгоритма. Суть алгоритма заключается в том, что в качестве информативных параметров ${\iota }_{n^{\text{'}}},$ ${\iota }_n$ выступают усеченные системы базисных ${\left\{\tilde{\beta }\right\}}_{n^{\text{'}}}\equiv {\left\{{\beta }_{\tilde{q}\text{'}}\right\}}_{n^{\text{'}}}$ и тестовых ${\left\{\tilde{\tau }\right\}}_n\equiv {\left\{{\tau }_{\tilde{q}}\right\}}_n$ функций $(\tilde{q}={\tilde{q}}_n\in \tilde{\mathcal{Q}}:$ $1\dots \tilde{Q}\mathrm{,}$ $\tilde{q}\text{'}={\tilde{q}}_{n^{\text{'}}}\in \tilde{\mathcal{Q}}\text{'}:$ $1\dots \tilde{Q}; \tilde{Q}\ll Q, \tilde{Q}\text{'}\ll Q^{\text{'}}).$ Ключом $\kappa =\kappa ({\iota }_n,{\iota }_{n^{\text{'}}})$ в данном случае выступают матрицы собственных и взаимных импедансов небольшой размерности ${\tilde{z}}_{n,n^{\text{'}}}=<{\tau }_{n,\tilde{q}},{\Psi }_{n,n^{\text{'}}}\left({\beta }_{n^{\text{'}},{\tilde{q}}^{\text{'}}}\right)>.$ Для рассматриваемой структуры усеченные СПФ целесообразно строить на основе следующих систем функций:

${\beta }_{n,q}\left(l\right)=\sqrt{2}\cos \left(\right(2q-1)\pi l/(2L\left)\right),{\tau }_{n,q}\left(l\right)={\beta }_{n,q}\left(l\right).$   (19)

Данные системы близки к системе СФ тонкого прямолинейного проводника длиной $L$ при небольших значениях отношений $L/\lambda$ [23]. Вычисления ${\tilde{z}}_{n,n^{\text{'}}}$ для рассматриваемой структуры целесообразно проводить в рамках обобщенного метода коллокаций с помощью выражений (9)–(11).

Алгоритм расчета состоит из следующих шагов:

• для текущих n, $n^{\text{'}}$ на основе информативных параметров ${\iota }_i, {\iota }_j$ соответствующих элементов вычисляется ключ $\kappa =\kappa ({\iota }_n,{\iota }_{n^{\text{'}}});$
• выполняется поиск ключа ${\kappa }^{\text{'}}$ в $\left\{\underset{¯}{\kappa }\right\},$ удовлетворяющего условию $\rho =|\kappa -{\kappa }^{\text{'}}|/\left|{\kappa }^{\text{'}}\right|\le {\rho }_{*},$ где ${\rho }_{*}$ – допустимое значение невязки ключей;
• если ключ ${\kappa }^{\text{'}}$ в k-й позиции вектора $\left\{\underset{¯}{\kappa }\right\},$ удовлетворяет условию $\rho \le {\rho }_{*}$ то $p_{n,n^{\text{'}}}^{\left(\overline{z}\right)}=k;$
• если ни один ключ ${\kappa }^{\text{'}}\in \left\{\underset{¯}{\kappa }\right\}$ не удовлетворяет условию $\rho \le {\rho }_{*},$ то вектор $\left\{\underset{¯}{\kappa }\right\}$ дополняется ключом $\kappa ,$ вектор $\left\{\underset{¯}{\overline{z}}\right\}$ – матрицей ${\overline{z}}_{n,n^{\text{'}}},$ а $p_{n,n^{\text{'}}}^{\left(\overline{z}\right)}=K^{\text{'}},$ где $K^{\text{'}}$ – число элементов дополненного вектора $\underset{¯}{\left\{\underset{¯}{\kappa }\right\}.}$

Алгоритм выполняется сначала для совпадающих, а затем – для несовпадающих $n$ и $n^{\text{'}}.$ На первом этапе определяется число матриц собственных импедансов, а на втором – число матриц взаимных импедансов в уникальном массиве матриц.

5. Исследование спектральных характеристик матрицы импедансов одиночного спирального элемента

Ранее было показано, что геометрия спирали (рис. 1, б) может быть определена с помощью шести параметров: радиусов нижнего и верхнего оснований $r_{\max }$ и $r_{\min },$ высоты спирали H, направления ее вращения s, числа витков $N_l$ и радиуса спирального проводника $\epsilon .$ При известных указанных параметрах по формулам (15) вычисляются коэффициенты $h_r$ и $h_z,$ входящие в уравнение образующей спирали (14). При исследовании спектральных характеристик матрицы импедансов одиночного спирального элемента целесообразным представляется использование нормированных параметров, где в качестве главного параметра, нормируемого к длине волны $\lambda$ следует использовать удвоенный максимальный радиус спирали $2r_{\max },$ определяющий максимальные габариты одиночного элемента при $H<2r_{\max }.$ Отношение $2r_{\max }/\lambda$ обозначим символом $\varsigma .$ При проведении численного моделирования исследовалась спираль, геометрия которой определялась следующими соотношениями: $2r_{\min }/\left(2r_{\max }\right)=$0,3682, $H/\left(2r_{\max }\right)=$0,025, $\epsilon /\left(2r_{\max }\right)=$0,0037, $N_l=\mathrm{2,}$$s=1.$ При этом $L/\left(2r_{\max }\right)\approx$4,323, или $L/\lambda =\chi \approx \mathrm{4,323}\varsigma .$ Последнее соотношение полезно знать в силу того, что в структурах, подобных рассматриваемым, возникают резонансы трех типов: тонкопроволочные, поверхностые и объемные. Тонкопроволочные резонансы определяются преимущественно длиной и формой проводника, являются высокодобротными, и наиболее ярко проявляются при полуцелых значениях $\chi ,$ не превышающих нескольких единиц.

При численном моделировании параметр $\varsigma$ менялся в пределах от 0,0125 до 0,53, при этом $\chi \in [\mathrm{0,054};\mathrm{2,28}].$ Матрицы импедансов формировались в рамках традиционного метода коллокаций, при этом число сегментов $M$ предполагалось равным 179, что при заданном отношении $\epsilon /\left(2r_{\max }\right)$ удовлетворяло условию (12). Здесь и далее волновое сопротивление среды $W_m$ полагалось равным $120\pi$ Ом, что справедливо для вакуума или воздушной среды. При этом предполагалось, что среда являлась недиссипативной, поэтому для вычисления волнового числа использовалось выражение $k=2\pi /\lambda .$

На рис. 3 представлены зависимости ${{\xi }^{\text{'}}}_n=\mathrm{Re}{\xi }_n,$ ${{{\xi }^{\text{'}}}^{\text{'}}}_n=\mathrm{sgn}\left(\mathrm{Im}\right({\xi }_n\left)\right)$ $\lg (1+|\mathrm{Im}\left({\xi }_n\right)\left|\right)$ (а) и $-\lg \left|{\xi }_n\right|$ от $\chi$ (б), индекс $n$ соответствует номеру кривой на графиках. Приведенные графики показывают, что зависимость СЧ от частоты имеет резонансный характер, и в рассматриваемом диапазоне мы наблюдаем четыре резонанса. Эти резонансы можно классифицировать как тонкопроволочные вследствие близкого соответствия значениям $n/2.$ При этом графики показывают, что резонанс для первой собственной функции наступает при $\chi ,$ несколько меньшем, чем $1/\mathrm{2,}$ а резонансы для последующих собственных функций возникают для значений $\chi ,$ больших, чем $n/\mathrm{2,}$ причем с ростом $n$ отклонения в большую сторону возрастают. Этот момент можно понять, если представить каждую собственную функцию в виде суперпозиции пары бегущих волн, распространяющихся по проводнику во взаимно противоположных направлениях. В этом случае становится очевидным, что фазовая скорость данных волн для $n=1$ будет меньше скорости света, а для $n>1$ – больше скорости света. Рост фазовой скорости, в свою очередь, связан с увеличением межвиткового взаимодействия, которому способствуют рост $\chi$ и особенности форм высших собственных функций. Также здесь можно отметить, что добротность резонансов довольно быстро падает с ростом n.

На рис. 4 приведен вид первых четырех нормированных собственных функций ${\stackrel{⌣}{\beta }}_n(l/L)/{\stackrel{⌣}{\beta }}_{n\mathrm{,}\max },$ вычисленных при соответствующих резонансных значениях $\chi .$ Здесь хорошо видно, что форма собственных функций близка к (19), но при этом общая амплитуда СФ при $n>1$ имеет существенную мнимую составляющую. Также здесь можно отметить асимметрию форм СФ, возрастающую с ростом n. Очевидно, что данная асимметрия определяется формой спирального элемента. Амплитудные распределения показывают, что СФ представляют собой преимущественно стоячую волну, при этом отличные от нуля значения тока в узлах при $n>1$ говорят о незначительном присутствии компоненты бегущей волны. Направление распространения данной компоненты на проводнике можно идентифицировать по направлению убывания амплитуды СФ.

 

Рис. 3. Зависимость ${{\mathbf{\xi }}^{\mathbf{\text{'}}}}_{\mathbf{n}}\mathbf{=}\mathit{R}\mathit{e}{\mathit{\xi }}_{\mathbf{n}}$ (сплошные кривые, а), ${{{\mathbf{\xi }}^{\mathbf{\text{'}}}}^{\mathbf{\text{'}}}}_{\mathbf{n}}\mathbf{=}\mathit{s}\mathit{g}\mathit{n}\mathbf{\left(}\mathit{I}\mathit{m}\mathbf{\right(}{\mathit{\xi }}_{\mathbf{n}}\mathbf{\left)}\mathbf{\right)}\mathbf{ }\mathbf{lg}\mathbf{(}\mathbf{1}\mathbf{+}\mathbf{|}\mathbf{Im}\mathbf{\left(}{\mathit{\xi }}_{\mathbf{n}}\mathbf{\right)}\mathbf{\left|}\mathbf{\right)}$ (штриховые кривые, а) и $\mathbf{-}\mathit{l}\mathit{g}\mathbf{\left|}{\mathit{\xi }}_{\mathbf{n}}\mathbf{\right|}$ (б) от $\mathit{\chi }\mathbf{;}$ номер маркера (линии) соответствуют номеру собственного значения

Fig. 3. Dependence ${{\mathbf{\xi }}^{\mathbf{\text{'}}}}_{\mathbf{n}}\mathbf{=}\mathit{R}\mathit{e}{\mathit{\xi }}_{\mathbf{n}}$ (solid curves, a), ${{{\mathbf{\xi }}^{\mathbf{\text{'}}}}^{\mathbf{\text{'}}}}_{\mathbf{n}}\mathbf{=}\mathit{s}\mathit{g}\mathit{n}\mathbf{\left(}\mathit{I}\mathit{m}\mathbf{\right(}{\mathit{\xi }}_{\mathbf{n}}\mathbf{\left)}\mathbf{\right)}\mathbf{ }\mathbf{lg}\mathbf{(}\mathbf{1}\mathbf{+}\mathbf{|}\mathbf{Im}\mathbf{\left(}{\mathit{\xi }}_{\mathbf{n}}\mathbf{\right)}\mathbf{\left|}\mathbf{\right)}$ (dashed curves, a) and $\mathbf{-}\mathit{l}\mathit{g}\mathbf{\left|}{\mathit{\xi }}_{\mathbf{n}}\mathbf{\right|}$ (b) from $\mathit{\chi }\mathbf{;}$ marker (line) number corresponds to number eigenvalue

 

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

6. Исследование эффективности алгоритма сжатия блочной матрицы и сходимости итерационного процесса

При численном моделировании структуры использовались геометрические параметры одиночного элемента, приведенные в предыдущем разделе. Считалось, что все элементы $L_{i_x\mathrm{,}{i}_y},$ связанные с одиночным элементом $L$ выражением (17), имеют один и тот же фиксированный угол ${\varphi }_{i_x,i_y}=\varphi ,$ который в расчетах полагался равным 0 градусов. Структура возбуждалась ПЭМВ W, ЭМП которой определяется выражениями (18). Амплитуда вектора электрического поля $E_0$ полагалась равной 1 В/м, начальная фаза ${\psi }_0$ – 0°, угол ${\varpi }_0$ – 0°, угол поляризации ${\varpi }_0=90$ (H-поляризация); угол падения ${\theta }_0$ принимал значения 0, 45 и 90°. Эти случаи мы будем называть случаями торцевого, бокового и нормального падения ПЭМВ соответственно.

 

Рис. 4. Вид первых четырех нормированных собственных функций ${\stackrel{\mathbf{⌣}}{\mathbf{\beta }}}_{\mathbf{n}}\mathbf{(}\mathit{l}\mathbf{/}\mathit{L}\mathbf{)}\mathbf{/}{\stackrel{\mathbf{⌣}}{\mathbf{\beta }}}_{\mathbf{n}\mathbf{,}\mathbf{m}\mathbf{a}\mathbf{x}}\mathbf{,}$ вычисленных при соответствующих резонансных значениях $\mathit{\chi }\mathbf{;}$ номера и составляющие собственных функций указаны на графиках

Fig. 4. Form of the first four normalized eigenfunctions ${\stackrel{\mathbf{⌣}}{\mathbf{\beta }}}_{\mathbf{n}}\mathbf{(}\mathit{l}\mathbf{/}\mathit{L}\mathbf{)}\mathbf{/}{\stackrel{\mathbf{⌣}}{\mathbf{\beta }}}_{\mathbf{n}\mathbf{,}\mathbf{max}}\mathbf{,}$ calculated for the corresponding resonant values $\mathit{\chi }\mathbf{;}$ numbers and components of eigenfunctions are indicated on the charts

 

Геометрия решетки определяется структуры четырьмя параметрами: расстояниями между соседними узлами решетки $h_x$ и $h_y$ вдоль соответствующих осей координат и числом элементов $N_x,$ $N_y$ вдоль соответствующих осей. Отношение $h_x/\lambda$ обозначим символом $\vartheta$. При численном моделировании предполагалось, что $h_x=h_y=h_a,$ $N_x=N_y=N_a.$ Рассматривались структуры с $N_a=\mathrm{2,4}\dots 24.$ Общее число элементов $N$ подобных структур определяется как $N_a^2.$ Для рассматриваемой структуры параметры $\varsigma$ и $\vartheta$ связаны соотношением: $\varsigma /\vartheta =$0,53, при этом данное соотношение, по сути, является коэффициентом заполнения структуры элементами. При $\varsigma /\vartheta =1$ соседние элементы структуры вплотную примыкают друг к другу. Исследования проводились в диапазоне $\vartheta \in [\mathrm{0,0625};1].$ Сходимость итерационного процесса в нашем случае зависит только от $N_a$ и $\vartheta .$ Коэффициент сжатия $K/N$зависит только от $N_a.$ На рис. 5 показан график зависимости коэффициента сжатия $K/N$ от числа элементов структуры N. Видно, что данная зависимость является линейной, и может быть приближенно описана функцией $K/N\approx$$(4/15)N=(4/15)N_a^2.$ Отметим, что при непосредственном расчете матричных элементов данная зависимость является квадратичной. Таким образом, предложенный алгоритм расчета элементов блочной матрицы позволяет существенно сократить затраты машинной памяти и машинного времени.

 

Рис. 5. Зависимость коэффициента сжатия $\mathit{K}\mathbf{/}\mathit{N}$ от числа элементов структуры $\mathit{N}$

Fig. 5. Dependence of the compression ratio $\mathit{K}\mathbf{/}\mathit{N}$ on the number of structure elements $\mathit{N}$

 

На рис. 6 представлена таблица сходимости итерационного процесса. Максимальное число итераций равно 100. Таким образом, достижение данного значения соответствует отсутствию сходимости итерационного процесса. С учетом указанной ранее связи между переменными $\sigma$ и $\chi ,$ тонкопроволочные резонансы на одиночных элементах структуры будут возникать при условии $\vartheta \approx \mathrm{0,217}n,$ в котором $n$ – номер резонанса. В результате для первых четырех собственных функций резонансы должны возникать в окрестности $\vartheta =$0,217, 0,434, 0,651, 0,868. При этом возможны смещения резонансных частот вследствие взаимодействия элементов структуры. Рисунок показывает, что ухудшение и отсутствие сходимости наблюдается в окрестности резонансных частот первых трех собственных функций, при этом первый резонанс не приводит к потере сходимости итерационного процесса, а резонанс четвертой собственной функции себя вообще никак не проявляет. Для более полного объяснения этого эффекта необходим расчет угловых распределений поля каждой СФ с определением пространственной ориентации максимумов поля.

 

Рис. 6. Таблица сходимости итерационного процесса: ось абсцисс – $\mathit{\vartheta }\mathbf{,}$ ось ординат – ${\mathit{N}}_{\mathbf{a}}\mathbf{,}$ значения в ячейках – число итераций

Fig. 6. Convergence table of the iterative process: x-axis – $\mathit{\vartheta }\mathbf{,}$ ordinate axis is ${\mathit{N}}_{\mathbf{a}}\mathbf{,}$ values in cells are number of iterations

 

7. Анализ и интерпретация результатов решения внутренней электродинамической задачи

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

${\eta }_n\left(l\right)=\sum _q{I}_q^{\left(n\right)}{\stackrel{⌣}{\beta }}_{n,q}\left(l\right),{\zeta }_n\left(l\right)=\sum _q{E}_q^{\left(n\right)}{\stackrel{⌣}{\beta }}_{n,q}\left(l\right);$

$l\in L,I_q^{\left(n\right)}\in i_n,E_q^{\left(n\right)}\in e_n.$

Для рассматриваемой структуры целесообразно оперировать векторами усредненных по $n$ амплитудных коэффициентов:

${\dot{f}}^{\left(a\right)}\supset {\dot{F}}_q^{\left(a\right)},{\dot{F}}_q^{\left(a\right)}=\frac{1}{N}\sum _n\left|F_q^{\left(n\right)}\right|,f\equiv i,e,F\equiv I,E.$

и соответствующими значениями невязок ${\rho }_n^{(f,a)}=|{\dot{f}}^{\left(a\right)}-f_n^{\left(a\right)}|/\left|{\dot{f}}^{\left(a\right)}\right|; f\equiv i,e.$ Расчеты показали, что в исследуемом частотном диапазоне для структур с $N_a=2\dots 24$ при углах падения ${\theta }_0=$0…90° средние значения этих невязок не превышают 0,09, а медианные значения меньше 0,025. Это, с одной стороны, подтверждает целесообразность рассмотрения усредненных векторов амплитудных коэффициентов, а с другой – говорит о высокой степени однородности токовых функций. Здесь и далее нормированные к максимальному значению векторы амплитудных коэффициентов будем обозначать как $f^{(a/\max )}=f^{\left(a\right)}/f_{\max }^{\left(a\right)}.$ Под $f_{\max }^{\left(a\right)}$ будем понимать максимальный амплитудный элемент вектора $f,$ а под $f_{nrm}^{\left(a\right)}=\left|f^{\left(a\right)}\right|$ – его евклидову норму. Супремумом вектора $f^{\left(a\right)}$ будем называть величину $f^{\left(\sup \right)}=\max (f_{\max }^{\left(a\right)},f_{nrm}^{\left(a\right)}).$ На рис. 7 сверху представлена нормированная спектрограмма ${\dot{e}}^{(a/\max )}\left(\vartheta \right)$ для случая $N_a=\mathrm{8,}$${\theta }_0=90°.$ Вдоль оси абсцисс спектрограммы отложены значения $\vartheta ,$ вдоль оси ординат – номера элементов вектора ${\dot{e}}^{(a/\max )}.$ Видно, что спектрограмма не зависит от $\vartheta ,$ при этом максимальная амплитуда наблюдается для четвертой собственной функции. Это обусловлено тем, что исследуемые спиральные элементы имеют число витков, равное двум. В целом можно отметить, что спектр является довольно широким в исследуемом диапазоне значений  На рис. 7 снизу представлена соответствующая нормированная спектрограмма ${\dot{i}}^{(a/\max )}\left(\vartheta \right)$ для аналогичного случая, а на рис. 8 – зависимости ${\dot{f}}_j^{(a/\sup )}={\dot{f}}_j^{\left(a\right)}\left(\vartheta \right)/{\dot{f}}^{\left(\sup \right)}$ $(f\equiv i,e, j\equiv \max ,nrm).$  Приведенные зависимости и спектрограмма показывают, что спектр собственных функций для рассматриваемой структуры в исследуемом диапазоне значений $\theta$ является узкополосным, и имеет при этом резонансный характер.

Разница уровней графиков ${\dot{f}}_{\max }^{(a/\sup )}$ и ${\dot{f}}_{nrm}^{(a/\sup )}$ $(f\equiv i,e)$ позволяет оценить степень вырожденности спектра. Так, если ${\dot{f}}_{\max }^{(a/\sup )}={\dot{f}}_{nrm}^{(a/\sup )},$ то спектр является вырожденным, и состоит только из одной спектральной компоненты. Таким образом, спектр стороннего поля является невырожденным во всем исследуемом диапазоне, а спектр токовых функций близок к вырожденному, и практически является таковым вблизи резонансных частот, определяющихся резонансными частотами для собственных значений элементов, из которых состоит метаструктура.

 

Рис. 7. Спектрограммы: сверху ${\dot{\mathbf{e}}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{m}\mathbf{a}\mathbf{x}\mathbf{)}}\mathbf{\left(}\mathit{\vartheta }\mathbf{\right)}\mathbf{,}$ – снизу – ${\dot{\mathbf{i}}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{m}\mathbf{a}\mathbf{x}\mathbf{)}}\mathbf{\left(}\mathit{\vartheta }\mathbf{\right)}\mathbf{;}$ ось абсцисс – $\mathit{\vartheta },$ ось ординат – индекс элемента векторов

Fig. 7. Spectrograms: top ${\dot{\mathbf{e}}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{m}\mathbf{a}\mathbf{x}\mathbf{)}}\mathbf{\left(}\mathit{\vartheta }\mathbf{\right)}\mathbf{,}$ – ${\dot{\mathbf{i}}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{m}\mathbf{a}\mathbf{x}\mathbf{)}}\mathbf{\left(}\mathit{\vartheta }\mathbf{\right)}\mathbf{;}$ below – $\mathit{\vartheta }\mathbf{,}$ abscissa axis –  ordinate axis – index of element of vectors

 

Рис. 8. Зависимости ${\dot{\mathbf{i}}}_{\mathbf{m}\mathbf{a}\mathbf{x}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{s}\mathbf{u}\mathbf{p}\mathbf{)}}$ (кривая 1), ${\dot{\mathbf{i}}}_{\mathbf{n}\mathbf{r}\mathbf{m}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{s}\mathbf{u}\mathbf{p}\mathbf{)}}$ (кривая 2), ${\dot{\mathbf{e}}}_{\mathbf{m}\mathbf{a}\mathbf{x}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{s}\mathbf{u}\mathbf{p}\mathbf{)}}$ (кривая 3) и ${\dot{\mathbf{e}}}_{\mathbf{n}\mathbf{r}\mathbf{m}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{s}\mathbf{u}\mathbf{p}\mathbf{)}}$ (кривая 4)

Fig. 8. Dependencies ${\dot{\mathbf{i}}}_{\mathbf{max}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{sup}\mathbf{)}}$ (curve 1), ${\dot{\mathbf{i}}}_{\mathbf{n}\mathbf{r}\mathbf{m}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{s}\mathbf{u}\mathbf{p}\mathbf{)}}$ (curve 2), ${\dot{\mathbf{e}}}_{\mathbf{m}\mathbf{a}\mathbf{x}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{s}\mathbf{u}\mathbf{p}\mathbf{)}}$ (curve 3) and ${\dot{\mathbf{e}}}_{\mathbf{nrm}}^{\mathbf{(}\mathbf{a}\mathbf{/}\mathbf{sup}\mathbf{)}}$ (curve 4) from

 

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

Заключение

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

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

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

Об авторах

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

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

Автор, ответственный за переписку.
Email: illuminator84@yandex.ru
ORCID iD: 0000-0002-9173-4936
SPIN-код: 9666-0814
ResearcherId: Q-9888-2017

доктор физико-математических наук, профессор кафедры физики 

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

Бассам Мохаммед-Али Аль-Нозайли

Самарский национальный исследовательский университет имени академика С.П. Королева

Email: bassam_91@mail.ru
SPIN-код: 7368-7223

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

Россия, 443086, Самара, Московское шоссе, 34

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

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