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

ВВЕДЕНИЕ

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

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

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

При разработке этих алгоритмов необходимо решить ряд электродинамических задач прикладной направленности.

Одна из таких задач – пересчет амплитудно-фазового распределения, измеренного на плоской поверхности сканирования P₀ (рис. 1), параллельной раскрыву рупорного облучателя, на плоскость входов фазовращателей проходной ФАР P₁ при различных углах поворота Θ облучающего рупора.

В более общем виде эту задачу можно сформулировать следующим образом: «пересчет известного амплитудно-фазового распределения поля на исходной плоскости P₀ на плоскость P₁, не параллельную исходной». Решение этой задачи позволяет избежать проведения дополнительных измерений поля на плоскости сканирования P^I₀, используя одно измерение, полученное при сканировании на плоскости P₀, для всех возможных углов отклонения Θ.

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

Для простоты будем рассматривать только линейную поляризацию вектора напряженности электрического поля и разместим его вдоль оси Оy, что не приведет к потере общности, так как результаты для случая линейной поляризации вектора E вдоль оси Ox и круговой поляризации могут быть получены поворотом системы координат (СК). Размеры плоскости сканирования P₀ должны быть таковы, чтобы при рассматриваемом положении облучателя амплитуда поля на ее краях составляла не более –30 дБ от максимального значения [1].

Пространственная модель для проведения сравнительного анализа численных методов приведена на рисунке 2.

Решетка диполей расположена в плоскости Oxy. P₀ – плоскость, параллельная решетке диполей, находящаяся на расстоянии z₀ = 5λ (где λ – длина волны) от начала координат. P₁ – плоскость, центр которой находится на расстоянии z₁ = 100λ от решетки диполей. Она повернута вокруг оси Оy на угол Θ_y, равный 10 градусам. Размеры плоскостей P₀ и P₁ одинаковы и равны 124,5λ на 127λ (250 на 255 точек с шагом 0,5λ).

Электрическое поле от модельной решетки диполей рассчитывается на плоскость P₀ (рис. 3). Затем с плоскости P₀ пересчитывается на повернутую плоскость P₁ каждым из исследуемых методов. Для проведения оценки результатов в качестве эталона используется поле, рассчитанное на повернутой плоскости P₁ непосредственно от решетки диполей. Размер модельной дипольной решетки 20 элементов (вдоль оси Ox) на 8 элементов (вдоль оси Oy) излучателей. Длина каждого диполя равняется половине длины волны. Расстояние между диполями – половина длины волны. Диполи ориентированы вдоль оси Oy, что соответствует вектору E вдоль оси Oy. Возбуждение решетки – равноамплитудное с перепадом фазы на 180 градусов при x = 0, что позволяет получить разностную диаграмму направленности (ДН) от поля на плоскости для оценивания точностных характеристик рассматриваемых численных методов.

1. Интерполяционный метод

Для расчета электрического поля на плоскости P₁, не параллельной исходной P₀ и повернутой к ней на угол Θ_y, можно воспользоваться интерполяционным методом (рис. 4).

Для реализации данного способа необходимо найти амплитудно-фазовое распределение (АФР) на каждой из промежуточных плоскостей P₁^Пр, P₂^Пр, ..., P_n^Пр, которое создается электрическим полем исходной плоскости P₀. Число плоскостей n зависит от расстояния между ними. Величина d_дис выбирается из условия, что расстояние между соседними точками пересечения искомой плоскости P₁ и параллельных плоскостей P₁^Пр, P₂^Пр, ..., P_n^Пр (узловыми точками) не должно быть больше чем λ/2.

Следовательно, d_дис ≤ λ sinΘy / 2. В общем случае, когда поворот плоскости P₁ происходит вокруг оси Oz или вокруг двух и более осей, необходимо произвести интерполяцию из узловых точек на регулярную сетку [2]. В рассматриваемом случае точки пересечения совпадают с регулярной сеткой плоскости P₁ и интерполяция не требуется.

2. Метод эквивалентных поверхностных токов (метод Гюйгенса – Кирхгофа)

Реализация данного метода основана на применении в расчетах электрического поля эквивалентных электрических и магнитных диполей. Электрические и магнитные поля заменяются соответственно на магнитные и электрические моменты.

Известны формулы [3] для электрического поля, создаваемого электрическим и магнитным диполем:

(1)

(2)

где эквивалентные электрический и магнитный моменты пропорциональны магнитному и электрическому полям соответственно [3]:

(3)

(4)

где H^S_m(x, y, z₀)вычисляется с использованием аппарата теории плоских волн [1], [4], [5].

В формулах (1), (2) диполи ориентированы вдоль оси Oz.

Так как рассматривается линейная поляризация вектора напряженности электрического поля E_y, то эквивалентный магнитный диполь также направлен вдоль оси Oy, а эквивалентный электрический диполь (соответствующий H_x) – вдоль оси Ox. Для учета этого условия в расчетах применяются повороты СК. Формулы (1)–(4) позволяют рассчитать поле на плоскости P₁, зная поле на плоскости P₀.

3. Метод поворотной трансформации спектра плоских волн

В данном методе задаются две СК [6]. Первая СК выступает в роли исходной СК (x, y, z), в которой находится излучающая плоскость P₀ с координатами (x, y, z₀). Во второй СК (xˆ, ŷ, zˆ) находится искомая плоскость P₁ с координатами (xˆ, ŷ, zˆ₁).

Единичные векторы r = (x, y, z) и rˆ = (xˆ, ŷ, zˆ) могут быть найдены с помощью матрицы поворота T следующим образом:

(5)

(6)

Матрица поворота T (7) вокруг трех осей в трехмерном пространстве представляет собой произведение трех матриц поворота вокруг каждой из осей, а в случае если поворот производится вокруг одной оси, то используется одна соответствующая повороту матрица: (8), (9) или (10).

(7)

где Wx(θx), Wy(θy), Wz(θz) – матрицы поворотов вокруг осей Ox, Oy и Oz соответственно, θ_x, θ_y, θ_z – углы поворота вокруг осей Ox, Oy и Oz против часовой стрелки для наблюдателя, смотрящего с конца оси.

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

(8)

(9)

(10)

При этом матрицы вращения обладают следующими характеристиками:

(11)

где W_ζ ^–1(χ) и W_ζ^T(χ)– обратная и транспонированная матрицы вращения соответственно.

В общем виде эти матрицы могут быть представлены следующим образом:

(12)

(13)

Если электрическое поле E(x, y, z) с волновым числом k задано на плоскости P₀, то электрическое поле на произвольных плоскостях, параллельных Oxy, задается в следующем виде [1], [6]:

(14)

где k_x, k_y, k_z – проекции волнового вектора, k_z = . а F(k_x, k_y) – спектр плоских волн электрического поля E(x, y, z), вычисляемый как:

(15)

Зададим новую СК так, чтобы ее оси были параллельны одноименным осям системы Oxy, расстояние от ее начала до плоскости P₁ равнялось нулю, а центр плоскости P₀ имел координаты (0, 0, –(z1 – z0)). Это позволит обнулить z в формуле (14), при этом в формуле (15) будет произведена подстановка d = –(z1 – z0).

Электрическое поле Eˆ (x, y, z) на повернутой плоскости P₁ находится путем подстановки (6) и (12) в (14):

(16)

Векторы f(kx, ky, kz) и fˆ(kˆx, kˆy, kˆz) в пространстве Фурье преобразуются один в другой также с помощью матриц T и T^–1. Поэтому, используя соотношение fˆ = Tf, аргумент экспоненты в выражении (16) можно упростить до 2πi(kˆ_x xˆ + kˆ_y ŷ).

При этом электрическое поле на повернутой плоскости P₁ будет определяться выражением:

(17)

где проекции волнового вектора k_x, k_y в уравнении (16) выражаются с использованием соотношения f = T^–1 fˆ следующим образом:

(18)

(19)

Замена в интеграле (16) переменных k_x, k_y на kˆ_x, kˆ_y осуществляется с помощью применения Якобиана:

(20)

где

(21)

С использованием (20) может быть реализовано преобразование электрического поля к виду:

(22)

где

(23)

Здесь

(24)

Графические иллюстрации по результатам моделирования процесса применения метода трансформации спектра приведены на рисунке 5.

4. Сравнение численных методов

Результаты анализа эффективности использования вышеописанных численных методов приведены на рисунке 6, где представлены азимутальные сечения ДН, рассчитанные с шагом 0,01 градуса.

Основное внимание уделено области одного из двух главных лепестков ДН и области равносигнального направления. Приведенные графики позволяют судить о степени приближения местоположения облучателя в пространстве к его оптимальному значению.

В таблице 1 приведены численные параметры полученных ДН.

Временные затраты на вычисления по каждому из методов были измерены при различном числе отсчетов на плоскостях P₀ и P₁. Соответствующий график представлен на рисунке 7. Из этого рисунка также видно, что метод интерполяции, как и ожидалось, предъявляет повышенные требования к объему оперативного запоминающего устройства (ОЗУ) ЭВМ. При достижении количеством элементов порога в 70 000 единиц происходит прекращение дальнейшего расчета из-за переполнения ОЗУ при вызове сильно расходующей память функции интерполяции. Объем доступной ОЗУ ЭВМ при выполнении численного моделирования составлял 11 ГБ.

ВЫВОДЫ

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