Solution of the internal and external problems of electrodynamics for a symmetrical vibrator with a reflector made of parallel straight conductors

Full Text

Введение

Моделирование распределения токов на многоэлементных излучающих структурах имеет особое значение в теории антенн. При проектировании таких излучающих структур необходимо учитывать взаимодействие между их отдельными элементами, что является довольно сложной задачей. Долгое время в инженерной практике для расчета взаимодействия использовался метод наведенных ЭДС. Основные принципы этого метода изложены в работе [1]. Метод наведенных ЭДС позволяет находить наведенные собственные сопротивления элементов антенны, а также амплитуды и фазы токов в пассивных элементах. Из недостатков метода отметим, что он накладывает определенные ограничения на длину и расстояние между элементами в силу использования приближенных распределений тока.

Сейчас для решения таких задач применяют системы автоматизированного проектирования, использующие в своей основе методы моментов [2], конечных элементов и конечных разностей [3]. К недостаткам такого подхода можно отнести высокие требования к ЭВМ, большую стоимость программного обеспечения, отсутствие математической модели в явном виде.

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

В [9] рассмотрен вопрос использования метода [5] для расчета распределений токов на директорной антенне [10]. Показано, что для директорной антенны в исследуемом диапазоне частот метод обладает хорошей сходимостью.

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

1. Внутренняя электродинамическая задача для многоэлементной структуры

Рассмотрим излучающую либо переизлучающую структуру, образованную совокупностью элементов $v_n$ $(n\in \mathcal{N}:1\dots N,$ $N$  – число элементов). Внутренняя задача электродинамики для рассматриваемой структуры сводится к операторной системе вида [9]:

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

В данной системе ${\Psi }_{n,n^{\text{'}}},$ $n\ne n^{\text{'}}$ – операторы взаимодействия элементов $v_n$ и $v_{n^{\text{'}}};$ ${\Psi }_{n,n}$ – собственный оператор элемента $v_n;$ ${\zeta }_n$ – функции, порождаемые сторонними полями, ${\eta }_{n^{\text{'}}}$ – токовые функции, подлежащие определению. К системе (1) применяется схема метода моментов [2]: с помощью систем базисных ${\left\{\beta \right\}}_{n^{\text{'}}}\equiv {\left\{{\beta }_{q^{\text{'}}}\right\}}_{n^{\text{'}}}$ и тестовых ${\left\{\tau \right\}}_n\equiv {\left\{{\tau }_q\right\}}_n$ функций $(q=q_n\in \mathcal{Q}:1\dots Q,$ $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{'}}}^{\left(n^{\text{'}}\right)}$ формируется система линейных алгебраических уравнений (СЛАУ) вида:

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

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

$\widehat{Z}\supset {\widehat{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)>,$ (4)

$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>.$

Здесь скобки $<>$ используются для обозначения скалярного произведения. Введем в рассмотрение массив уникальных матриц $\left\{{\underset{¯}{\widehat{z}}}_k\right\}\equiv \left\{\underset{¯}{\widehat{z}}\right\}$  $(k\in \mathcal{K}:1\dots K).$ Массив $\left\{p_{n,n^{\text{'}}}^{\left(\widehat{z}\right)}\right\}$ с элементами $p_{n,n^{\text{'}}}^{\left(\widehat{z}\right)}=k$ будет устанавливать соответствие между парой индексов $n,n^{\text{'}}\in \mathcal{N}$ и индексом $k\in \mathcal{K}.$ Таким образом, имеем:

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

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

Более эффективные методы расчета можно построить, опираясь на информативные параметры ${\iota }_n$ и ${\iota }_{n^{\text{'}}}$ n-го и $n^{\text{'}}$-го элементов структуры. Здесь нужно ввести понятие ключа $\kappa ,$ вычисляемого на основе информативных параметров: $\kappa =\kappa ({\iota }_n,{\iota }_{n^{\text{'}}}).$ Ключ $\kappa$ должен однозначно определять пары взаимодействий. Вектор уникальных ключей обозначим как $\left\{\underset{¯}{\kappa }\right\}.$ Процедура расчета массивов $\left\{p^{\left(\widehat{z}\right)}\right\}$ и $\left\{\underset{¯}{\widehat{z}}\right\}$ в данном случае состоит из следующих шагов:

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

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

2. Прямые и итерационные методы решения СЛАУ

Заключительным этапом решения внутренней задачи является решение блочной СЛАУ (3). Решать СЛАУ прямыми методами (методы LU-разложения, вращения, отражения, Гаусса и т. д. [6]) в случае использования базисных и тестовых функций универсального вида может быть нецелесообразным в силу ряда причин. Во-первых, даже для структур, содержащих не более сотни элементов, размерность матрицы будет довольно большой, а сложность решения пропорциональна третьей степени размерности. Также в данном случае будет наблюдаться накопление погрешностей округления [6]. Во-вторых, приведенный ранее алгоритм формирования внутренней задачи предполагает использование методов решения СЛАУ, не приводящих к модификации исходной блочной матрицы в процессе расчета или к выделению дополнительной памяти в размере, необходимом для хранения полной исходной блочной матрицы.

В качестве одного из альтернативных вариантов можно использовать итерационные методы, основанные на расщеплении матрицы (методы типа последовательной верхней релаксации [6]), однако в случае отсутствия диагонального преобладания в матрице $\widehat{Z}$ они не обеспечивают сходимости итерационного процесса.

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

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

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

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

${\widehat{y}}_n={\widehat{z}}_{n,n}^{-1},{\widehat{w}}_{n,n^{\text{'}}}={\widehat{p}}_n{\widehat{z}}_{n,n^{\text{'}}}$ (7)

– обращенные собственные матрицы и весовые матрицы соответственно;  выступают в качестве матриц предобуславливателя.

В процедуре (6) можно также использовать более простой вариант

${\widehat{p}}_n=diag{\left({\widehat{z}}_{n,n^{\text{'}}}\right)}^{-1}.$ (8)

Здесь оператор $diag\left(\widehat{z}\right)$ предполагает построение диагональной матрицы, размер которой совпадает с $\widehat{z},$ а в качестве диагональных элементов используются элементы $z_{q,q}\in \widehat{z}.$ Если положить во второй сумме (6) $l+1=l,$ то формула будет соответствовать методу простой итерации. Критерий оценки сходимости строится в соответствии с неравенством

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

где $\xi$ – сколь угодно малое наперед заданное число.

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

Многоэлементная тонкопроволочная структура $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\left(l\right),$ $Y_n\left(l\right),$ $Z_n\left(l\right)$ – гладкие функции. Под $L_n=$ $L_{n,\max }-L_{n,\min }$также будем понимать длину n-го проводника. Интегральное представление электромагнитного поля (ИП ЭМП) такой структуры можно записать в виде [12]:

$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;$ (10)

здесь $I_{n^{\text{'}}}\left(l^{\text{'}}\right)$ – распределение полного тока по образующей $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{'}}}\left(l^{\text{'}}\right)$ целесообразно представление в виде рядов типа (2). Исходное ИП ЭМП (10) при этом приобретет вид

$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{'}},$ (11)

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

${\mathbf{E}}^{\text{(ext)}}\left({\mathbf{r}}_{\mathrm{n}}\right(l\left)\right)+\mathbf{E}\left({\mathbf{r}}_{\mathrm{n}}\right(l\left)\right))\cdot {\hat{\mathbf{l}}}_{\mathrm{n}}(l)=0.$ (12)

Умножая поочередно (12) на тестовые функции ${\tau }_{n,q}\left(l\right)$ и интегрируя по l, получаем СЛАУ для вычисления $I_{q^{\text{'}}}^{\left(n^{\text{'}}\right)},$ по форме совпадающую с (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}$ (13)

$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),$ (14)

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

Конкретизируем способ вычисления интегралов в (13) с помощью процедуры сегментации проводников, предполагающей представление n-го проводника в виде совокупности $M+1$ узлов $L_n^{\left(M\right)}:$ $r_{n\mathrm{,1}},$ $r_{n\mathrm{,2}},$ $\dots ,$  $r_{n,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}),$ (15)

где ${\dot{l}}_{n,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}).$ (16)

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

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

$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{'}})},$ (17)

$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{'}},$ (18)

${\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}},$ (19)

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

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

4. Анализируемая излучающая структура

Рассматриваемая в статье структура является одним из вариантов панельной антенны, представляющей собой активный вибратор $L_a\equiv L_1,$ расположенный на высоте $z_{\max }$ над решеткой однотипных пассивных вибраторов (рефлекторов) $L_{r,r}\equiv L_{r+1}$ $(r=1\dots N_r).$ Геометрия структуры приведена на рис. 1. Все вибраторы являются тонкими и идеально проводящими, имеют одинаковый диаметр $2\epsilon \ll \lambda$ и расположены симметрично относительно плоскости zOx. Также исследуемая структура обладает симметрией относительно плоскости yOz. Длину активного вибратора далее будем обозначать как $L_a,$ длину рефлекторов – как $L_r.$ Распределения токов обозначим в соответствии с обозначениями самих вибраторов как $I_n\left(l\right)$ $(n\in \mathcal{N}:1\dots N=$ $N_r+1).$ Размер рефлектора вдоль оси $Ox$ равен $2x_{\max },$ таким образом, общий размер рефлектора равен $2x_{\max }\times L_r.$

В центре активного вибратора $L_a$ имеется разрыв длиной $b\ll \lambda ,$ в который помещен генератор сторонней ЭДС g. Под действием генератора сторонней ЭДС, имеющего амплитуду, равную $U$ В, на активном вибраторе возникает касательная компонента стороннего электрического поля $E^{\left(ext\right)}\left(l\right),$ равная нулю всюду за исключением области зазора, где она равна $U/b.$ Эта компонента создает на активном вибраторе распределение тока $I_1\left(l\right),$ который, в свою очередь, создает токи на пассивных элементах антенны.

Образующая любого проводника структуры может быть описана выражением:

$L_n:r\left(l\right)=\widehat{x}{x}_n+\widehat{y}l+\widehat{z}{z}_n,l\in [-L_n/2;L_n/2],$ (21)

здесь $x_n,$  $z_n$ – координаты точки перечечения n-го проводника с плоскостью xOz; для активного вибратора   $x_1=\mathrm{0,}$ $z_1=z_{\max };$ для проводников рефлектора $x_r=2x_{\max }(r-1)/$ $(N_r-1)-x_{\max },$ $z_r=\mathrm{0,}$  $r$ – номер проводников рефлектора.

 

Рис. 1. Геометрия исследуемой структуры

Fig. 1. Geometry of the structure under study

 

Используем ключевой алгоритм расчета массивов $\left\{\underset{¯}{\widehat{z}}\right\}$ и $\left\{p^{\left(\widehat{z}\right)}\right\},$ описанный в первом разделе. В данном случае нам необходимо дополнить его правилами выбора информативных параметров ${\iota }_n,$ расчета ключей и сравнения ключей $\kappa .$ Рассматриваемая структура имеет плоскость $xOz$ в качестве плоскости симметрии, а входящие в нее элементы создают азимутально-независимые электромагнитные поля, могут отличаться только длиной $(L_a$ или $L_r)$ и координатами $x_n,$ $z_n.$ Таким образом, информативные параметры n-го элемента представляют собой массив ${\iota }_n=\{L_n,x_n,z_n\}.$ Ключ $\kappa$ для элементов $n$ и $n^{\text{'}}$ может быть вычислен следующим образом:

$\kappa =\{L_n,L_{n^{\text{'}}},d\},d=\sqrt{{(x_n-x_{n^{\text{'}}})}^2+{(z_n-z_{n^{\text{'}}})}^2}.$ (22)

После сегментации такой структуры к ней в полной мере становятся применимы выражения (3) и (17). Длины проводников подберем так, чтобы они состояли из равных сегментов длиной ${\Delta }_a=L_a/M_a,$ ${\Delta }_r=L_r/M_r,$ здесь $M_a$ и $M_r$ – число сегментов активного вибратора и число сегментов рефлекторов соответственно. Длина сегмента выбирается в соответствии с (20). С учетом симметрии структуры относительно плоскости $zOx$ целесообразным является четное число сегментов на элементе: $M=2\tilde{M}.$ В этом случае конец сегмента с индексом $\tilde{M}$ будет началом сегмента с индексом $\tilde{M}+1$, а соответствующая точка будет принадлежать плоскости симметрии. Ширину зазора b, в который помещается генератор g, после сегментации можно принять равной $2\Delta .$ Далее при моделировании будем использовать СПФ, построенные на основе (18) и соответствующие различным вариантам метода коллокаций.

1. СПФ первого типа, учитывающие симметрию и способ возбуждения структуры:

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

${\tau }_{n,q}\left(l\right)={\beta }_{n,q}\left(l\right),q\in \mathcal{Q}:1\dots Q=\tilde{M};$

Выражения (17) при этом приобретают вид:

$Z_{q,q^{\text{'}}}^{(n,n^{\text{'}})}={\mu }_{q,q^{\text{'}}}^{(n,n^{\text{'}})}+{\mu }_{q,\tilde{q}\text{'}}^{(n,n^{\text{'}})}+{\mu }_{\tilde{q},q^{\text{'}}}^{(n,n^{\text{'}})}+{\mu }_{\tilde{q},\tilde{q}\text{'}}^{(n,n^{\text{'}})},$ (24)

$E_q^{\left(n\right)}={\nu }_q^{\left(n\right)}+{\nu }_{\tilde{q}}^{\left(n\right)}.$

В выражениях (23) и (24) $\tilde{q}=2\tilde{M}+1-q,$ ${\tilde{q}}^{\text{'}}=$$2{\tilde{M}}^{\text{'}}+1-q^{\text{'}}.$

2. СПФ второго типа, учитывающие симметрию и способ возбуждения структуры, а также вид собственных функций на ее отдельных элементах [11]:

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

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

При проведении численного моделирования в качестве главного параметра, нормируемого к длине волны, была выбрана длина активного вибратора $L_a,$ отношение $L_a/\lambda$ в дальнейшем будем обозначать как $\gamma .$ Геометрия исследуемой структуры определялась следующими нормированными соотношениями: $z_{\max }/L_a=$ $\mathrm{0,78125,}$ $L_r/L_a=$ $2x_{\max }/L_a=\mathrm{3,125,}$ $\epsilon /L_a=$$\mathrm{1,875}\cdot {10}^{-3}.$ Половинное число сегментов на активном вибраторе и элементах рефлектора предполагалось равным ${\tilde{M}}_a=116$ ${\tilde{M}}_r=119$ соответственно. Число проекционных функций для второго варианта СПФ было равно 12 на активном вибраторе и 5 на элементах рефлектора. Амплитуда  генератора ЭДС полагалась равной 1 В. Моделирование проводилось в диапазоне $\gamma =\mathrm{0,25}\dots \mathrm{0,55.}$ Отметим, что нормированная ширина зазора $d/L_a$ оказалась равной $1/116.$ При исследовании сходимости итерационных процессов в (9) предполагалось, что $\xi ={10}^{-3}.$ Исследования проводились для различного числа элементов рефлектора $N_r$ – от 5 до 19, которые предполагались нечетными. При $N_r=19$ в случае СПФ первого типа матрица $\widehat{Z}$ имела размер $2377\times 2377$ элементов, в случае СПФ второго типа – $107\times 107$ элементов. Очевидно, что в последнем случае решение СЛАУ прямыми методами не представляет особых проблем с точки зрения затрат машинной памяти и машинного времени.

На рис. 2 представлена визуализация структуры блочных матриц для различного числа элементов рефлектора. В центре каждого блока обозначено значение $p_{n,n^{\text{'}}}^{\left(\widehat{z}\right)},$ оси абсцисс представляют собой индексы столбцов матрицы $\widehat{Z},$ оси ординат – индексы строк матрицы $\widehat{Z}.$ На цветовой шкале отложены индексы вектора $\left\{\underset{¯}{\kappa }\right\}.$ Из представленного рисунка видно, что число элементов вектора $\left\{\underset{¯}{\kappa }\right\}$ существенно меньше, чем $N^2.$ Анализ для рефлекторов с нечетным числом элементов показал, что $K=2N,$ то есть сложность расчета блочной матрицы с учетом представленного алгоритма, имеет оценку $O\left(N\right),$ а не $O\left(N^2\right),$ что справедливо в случае непосредственного расчета блочной матрицы. Таким образом, предложенный алгоритм расчета матричных элементов демонстрирует свою эффективность. В частных случаях можно практически в два раза уменьшить K. Так, в случае, когда у элементов рефлектора и активного вибратора сегменты равной длины, справедливо тождество ${\widehat{z}}_{n,n^{\text{'}}}={\widehat{z}}_{n^{\text{'}},n}^T,$ которое можно использовать при формировании алгоритма расчета элементов матрицы $\widehat{Z}.$

 

Рис. 2. Структура блочных матриц в случае 5 (а), 9 (б), 15 (в) и 19 (г) элементов рефлектора; в центре блока обозначен номер вектора существенных элементов, оси абсцисс – индексы столбцов матриц, оси ординат – индексы строк матриц; на цветовой шкале отложены номер вектора существенных элементов

Fig. 2. Structure of block matrices in the case of 5 (a), 9 (b), 15 (c) and 19 (d) reflector elements; in the center of the block the number of the vector of essential elements is indicated, the abscissa axes are the indices of the matrix columns, the ordinate axes are the indices of the matrix rows; the number of the vector of essential elements is shown on the color scale

 

На рис. 3 показаны результаты исследования зависимости сходимости итерационного процесса (6) от $\gamma$ при различном числе элементов рефлектора (указаны в обозначениях линий) для СПФ первого (а) и второго (б) типа. Как можно видеть из представленных графиков, итерационный процесс расходится для обеих СПФ при $N_r>9$ в окрестности $\gamma \approx \mathrm{0,47,}$ отсутствие сходимости также наблюдается в окрестности $\gamma \approx \mathrm{0,3}$ при $N_r>15$ для СПФ первого типа и при $N_r>17$ для СПФ второго типа. Отметим, что $\gamma \approx \mathrm{0,3}$ соответствует второму резонансу одиночного проводника рефлектора, а $\gamma \approx \mathrm{0,47}$ – первому резонансу активного вибратора. Также здесь можно видеть, что при $N_r\ge 13$ происходит ухудшение сходимости итерационного процесса при малых значениях $\gamma ,$ связанное с приближением к точке первого резонанса проводников рефлектора $(\gamma \approx \mathrm{0,15}).$ Для директорной антенны, рассмотренной в [9], сходимость итерационного процесса наблюдалась при любых значениях $\gamma .$ Таким образом, можно сделать вывод, что отсутствие сходимости итерационного процесса для рассматриваемой структуры возникает при очень близком расположении элементов рефлектора, приводящем к усилению межэлементного взаимодействия. Резонансные явления в проводниках лишь усугубляют эту проблему. В нерезонансных случаях итерационный подход остается эффективным и позволяет за 3–30 итераций получить решение требуемой точности, что в вычислительном плане существенно проще прямого решения СЛАУ, которое можно использовать в резонансных случаях для гарантированного получения решения. В дальнейшем решение проблемы сходимости для рассматриваемой структуры авторы связывают с введением в итерационный алгоритм предобуславливающей матрицы, учитывающей резонансные эффекты.

 

Рис. 3. Исследование сходимости итерационного процесса решения внутренней задачи для СПФ первого (а) и второго (б) типа; номер кривой соответствует числу элементов рефлектора

Fig. 3. Study of the convergence of the iterative process of solving the internal problem for SPF of the first (a) and second (b) types; the number of the curve corresponds to the number of elements

 

На рис. 4 представлено сравнение амплитудных распределений тока на активном вибраторе, полученных при $N_r=9$ с помощью прямой и итерационной процедуры с использованием СПФ первого и второго типов. Из графиков можно видеть, что распределения тока имеют лишь незначительные визуальные отличия, что подтверждает корректность выбора числа проекционных функций в рамках СПФ второго типа.

 

Рис. 4. Сравнение амплитудных распределений тока на правом плече активного вибратора, полученных различными методами для $\gamma =\mathrm{0,24}$  (а) и $\gamma =\mathrm{0,48}$  (б): 1 – строгий метод, СПФ первого типа; 2 – итерационный метод, СПФ первого типа; 3 – итерационный метод, СПФ второго типа; $N_r=9$

Fig. 4. Comparison of amplitude current distributions on the right arm of the active vibrator, obtained by various methods for $\gamma =\mathrm{0,24}$  (a) and $\gamma =\mathrm{0,48}$  (b): 1 – strict method, SPF of the first type; 2 – iterative method, SPF of the first type; 3 – iterative method, SPF second type; $N_r=9$

 

На рис. 5 представлены зависимости нормированных максимальных амплитуд токов $A=$$\lg (I_{\max ,j}/I_{\max ,a})$ от $\gamma$ на элементах структуры для случаев $N_r=9$ и $N_r=19.$ Индекс $a$ на графиках и в формуле для $A$ соответствует активному вибратору, числовой индекс – номеру элемента рефлектора. На рисунках приведены кривые не для всех элементов рефлектора в силу симметрии структуры относительно плоскости yOz. В случае $N_r=9$ под активным вибратором расположен пятый элемент рефлектора, а в случае $N_r=19$ – десятый элемент. Как видно из представленных графиков, резонансы элементов рефлектора никак себя не проявляют, в отличие от первого резонанса на активном вибраторе $(\gamma =\mathrm{0,47}).$ Амплитуды токов на рефлекторах меньше амплитуды тока на активном вибраторе более чем на порядок.

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

 

Рис. 5. Зависимость нормированных максимальных амплитуд токов на активном вибраторе (кривая «a») и элементах рефлектора (кривые 1–5, 1–10) от $\gamma$  для рефлектора, содержащего 9 (а) и 19 (б) элементов

Fig. 5. Dependence of the normalized maximum amplitude of currents on the active vibrator (curve «a») and reflector elements (curves 1–5, 1–10) on $\gamma$  for a reflector containing 9 (a) and 19 (b) elements

 

На рис. 6, а приведен анализ зависимости корреляции форм амплитудных распределений тока на соседних элементах рефлектора от $\gamma ,$ выполненный по формуле

${\zeta }_r=\frac{{\displaystyle {\int }_{L_r}\left|I_r\right(l\left)\right|\left|I_{r+1}\right(l\left)\right|dl}}{\sqrt{{\displaystyle {\int }_{L_r}\left|I_r\right(l){|}^{2}dl}}\sqrt{{\displaystyle {\int }_{L_r}\left|I_{r+1}\right(l){|}^{2}dl}}},r=\mathrm{1...}{N}_r,$ (26)

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

${\rho }_r=2\frac{{\displaystyle {\int }_{L_r}\left(\left|\right|I_r\left(l\right)|-|I_{r+1}\left(l\right)\left|\right|\right)dl}}{{\displaystyle {\int }_{L_r}\left(\left|I_r\right(l\left)\right|+\left|I_{r+1}\right(l\left)\right|\right)dl}},r=\mathrm{1...}{N}_r.$ (27)

 

Рис. 6. Зависимость степени корреляции ${\zeta }_r$  (а) и относительной невязки амплитуд ${\rho }_r$  (б) на соседних проводниках рефлектора от $\gamma ;$  рефлектор содержит 9 элементов, значение  соответствует номеру кривой

Fig. 6. Dependence of the degree of correlation ${\zeta }_r$  (a) and the relative discrepancy of amplitudes ${\rho }_r$  (b) on adjacent conductors of the reflector on $\gamma ;$  reflector contains 9 elements, the value corresponds to the curve number

 

Соответствующие зависимости показаны на рис. 6, б. Здесь видно, что в нерезонансном случае наименьшую разницу амплитуд имеют соседние элементы рефлектора, расположенные под активным вибратором, при этом в резонансном случае рассматриваемая разность существенно увеличивается. В целом по результатам анализа графиков, представленных на рис. 5 и рис. 6, можно сделать вывод о гладкости рассмотренных зависимостей от $\gamma$ в нерезонансных случаях, которая нарушается в случае приближения к резонансу.

На рис. 7 приведена зависимость входного сопротивления активного вибратора от $\gamma ,$ рассчитанная при различных значениях $N_r.$ Графики показывают, что наличие рефлектора в данном случае практически не влияет на входное сопротивление, которое практически полностью совпадает с аналогичными зависимостями для одиночного тонкого симметричного вибратора. Резонанс достигается при $\gamma \approx \mathrm{0,47,}$ при этом значение входного сопротивления становится чисто активным и равным примерно 73 Ом.

 

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

Fig. 7. Dependence of the real and imaginary parts of the input resistance of the active vibrator for different numbers of reflector elements, the number of elements and the input resistance component are indicated in the designation of the curves

 

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

На рис. 8 представлены результаты сравнения нормированных меридианных диаграмм направленности (ДН), полученных различными методами для $\gamma =\mathrm{0,24}$ (а), $\gamma =\mathrm{0,32}$ (б), и $\gamma =\mathrm{0,48}$ (в). Число элементов рефлектора было равно девяти. Во-первых, здесь можно увидеть, что результаты, полученные различными методами, практически не имеют визуальных отличий. Во-вторых, при $\gamma <\mathrm{0,32}$ рефлектор выполняет свою функцию, блокируя значительную часть излучения в нижнее полупространство. При $\gamma =\mathrm{0,48}$ данное свойство рефлектора теряется и происходит формирование лепестка обратного излучения.

На рис. 9 представлены результаты расчета нормированных меридианных ДН при тех же значениях $\gamma$ для структур с рефлекторами, имеющими различное число элементов. Здесь видно, что увеличение числа элементов рефлектора при используемых параметрах моделирования ведет к существенному снижению лепестка обратного излучения, при этом в случае $\gamma =\mathrm{0,48}$ вместе с уменьшением уровня обратного лепестка происходит его сужение. Также здесь можно заметить, что с ростом $\gamma$ скорость падения амплитуды обратного лепестка при увеличении числа элементов рефлектора уменьшается.

 

Рис. 8. Сравнение нормированных амплитудных меридианных ДН $(\phi =0),$  полученных различными методами для $\gamma =\mathrm{0,24}$  (а), $\gamma =\mathrm{0,32}$  (б), и $\gamma =\mathrm{0,48}$  (в): 1 – строгий метод, СПФ первого типа; 2 – итерационный метод, СПФ первого типа; 3 – итерационный метод, СПФ второго типа; $N_r=9$

Fig. 8. Comparison of normalized amplitude meridian patterns $(\phi =0),$  obtained by various methods for $\gamma =\mathrm{0,24}$  (a), $\gamma =\mathrm{0,32}$  (b), and $\gamma =\mathrm{0,48}$  (c): 1 – strict method, SPF of the first type; 2 – iterative method, SPF of the first type; 3 – iterative method, SPF of the second type; $N_r=9$

 

Рис. 9. Сравнение нормированных амплитудных меридианных ДН $(\phi =0)$  для $\gamma =\mathrm{0,24}$  (а), $\gamma =\mathrm{0,32}$  (б) и $\gamma =\mathrm{0,48}$  (в), рассчитанных для различного числа элементов рефлектора (указаны в обозначении кривых)

Fig. 9. Comparison of normalized amplitude meridian patterns $(\phi =0)$  for $\gamma =\mathrm{0,24}$  (a), $\gamma =\mathrm{0,32}$ (b), and $\gamma =\mathrm{0,48}$  (c), calculated for different numbers of reflector elements (indicated in the captions to the curves)

 

Заключение

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

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

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

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

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

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

Анализ зависимости входного сопротивления активного вибратора от  показал, что наличие рефлектора в данном случае практически не влияет на входное сопротивление, которое практически полностью совпадает с аналогичными зависимостями для одиночного тонкого симметричного вибратора.

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

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

About the authors

Dmitry P. Tabakov

Povolzhskiy State University of Telecommunications and Informatics

Author for correspondence.
Email: illuminator84@yandex.ru
ORCID iD: 0000-0002-9173-4936
SPIN-code: 9666-0814
ResearcherId: Q-9888-2017

Doctor of Physical and Mathematical Sciences, professor of the Department of Radioelectronic Systems

Russian Federation, 23, L. Tolstoy Street, Samara, 443010

Bassam Mohammed-Ali Al-Nozaili

Samara National Research University

Email: bassam_91@mail.ru
SPIN-code: 7368-7223

graduate student, assistant professor of the Department of Physics

Russian Federation, 34, Moskovskoye shosse, Samara, 443086

References

Supplementary files