CC BY
66
А.Г. Дмитренко, В.А. Колчин
ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ РАССЕЯНИЯ ЭЛЕКТРОМАГНИТНЫХ ВОЛН НА ТРЕХМЕРНЫХ ИДЕАЛЬНО ПРОВОДЯЩИХ ТЕЛАХ
Предложен и реализован вариант метода дискретных источников для численного решения задач электромагнитного рассеяния на структурах, составленных из конечного числа трехмерных идеально проводящих тел, ограниченных гладкими поверхностями произвольной формы. Кратко описаны возможности разработанного пакета прикладных программ. Приведены некоторые результаты численных расчетов, касающиеся влияния электромагнитного взаимодействия тел на величины бистатических сечений рассеяния образованных ими структур.
Значительный интерес для исследователей представляет изучение рассеяния радиоволн совокупностью идеально проводящих тел, имеющих размеры, сравнимые с длиной волны. Этот интерес обусловлен необходимостью решения ряда практически важных проблем, например таких, как проблемы электромагнитной совместимости, конструирования многоэлементных антенных систем и др.
Корректная постановка исследований подобного рода приводит к необходимости решения граничных задач теории рассеяния на системах трехмерных взаимодействующих (в электромагнитном смысле) тел. Известно относительно небольшое число работ, посвященных решению таких задач. В редких случаях [1] удается получить аналитическое решение. В большинстве же случаев приходится прибегать к численным методам. Последние могут быть основаны как на использовании уравнений Максвелла в дифференциальной форме [2], так и на использовании интегральных соотношений теории электромагнитного поля [3]. Однако соответствующие вычислительные алгоритмы получаются чрезвычайно емкими по затратам ресурсов и времени ЭВМ. Для конечных методов это обусловлено необходимостью распространения вычислений на всю рассматриваемую область пространства, для методов интегральных уравнений - необходимостью вычисления большого количества поверхностных или объемных интегралов.
В последние годы, применительно к решению задач электромагнитного рассеяния, получил существенное развитие метод дискретных источников [4]. В этом методе неизвестное поле в рассматриваемой области и на ее границах представляется в виде конечной линейной комбинации полей некоторой системы источников, размещенных вне рассматриваемой области. Такая конструкция удовлетворяет системе уравнений Максвелла и условиям излучения (где это необходимо). Коэффициенты линейной комбинации определяются путем удовлетворения граничным условиям на поверхности рассеивателей. В силу своей идейной простоты метод удобен в качестве основы для построения решений задач электромагнитного рассеяния на телах самой общей формы, не обладающих симметрией вращения. В [5 - 10] развиты варианты метода дискретных источников для решения задач рассеяния на одиночных неосесимметричных телах с различными физическими свойствами. Накопленный опыт позволяет сделать соответствующие обобщения на системы взаимодействующих тел.
В данной работе численный метод решения задач электромагнитного рассеяния на одиночном идеально проводящем теле, предложенный в [5, 6], обобщен на случай системы взаимодействующих идеально проводящих тел, каждое из которых является неосесимметричным и имеет форму, отличную от формы других тел, входящих в систему. Даны математическая формулировка метода и краткое описание возможностей реализованной на его основе программы для расчета компонент рассеянного поля. Приведены некоторые результаты расчетов методического характера, а также результаты исследований влияния взаимодействия между телами на величины бистатических сечений рассеяния образованных ими структур.
ФОРМУЛИРОВКА ЗАДАЧИ И МЕТОД ЕЕ РЕШЕНИЯ
Геометрия задачи показана на рис. 1. Будем рассматривать стационарную (зависимость от времени выбрана в виде ехр(-/'юґ)) задачу дифракции электромагнитного поля {¿0, Н0} на структуре, состоящей из
Q непересекающихся идеально проводящих тел Вц, ограниченных поверхностями Бде С(2, а) (д = 1,2,..., Q),
3 Q
размещенной в однородной среде Бе = Я \( и О ) с
д=1
диэлектрической и магнитной проницаемостями єе и це в декартовой системе координат с центром, выбранным, для определенности, внутри рассеивателя с порядковым номером д = 1. Требуется найти рассеянное поле {Ёе, Не} в области Бе.
Математическая постановка задачи имеет следующий вид:
УхЁе = т\хеИе, УхИе =-1Ю£еЁе в Бе; (1)
Пд х Ее = -пд х Е0 на Бд , q = 1,0; (2)
{&Ёе ;Тм? е }х Я / Я +
+ {4^еИе;-фГеЁе} = 0(Я~'), Я , (3)
где Пд - единичные векторы нормалей к поверхностям
о о о 1 /о
Бд, Я = (х + у + 2) , а х Ь - векторное произведение.
Суть предлагаемого метода заключается в следующем. Введем (рис. 1) внутри каждого из рассеивателей Вд вспомогательную поверхность Бед=КедБд, подобную поверхности рассеивателя в смысле гомотетии с центром в точке 0д. Если поверхность Бд является центральной, центр гомотетии выбирается совпадающим с центром поверхности. Коэффициенты
гомотетии (подобия) Ке, д (д = 1,е) характеризуют удаление вспомогательных поверхностей от поверхностей соответствующих тел, их значения лежат в интервале 0<Ке, д<1 (при Ке, д=0 вспомогательная поверхность стягивается в точку, при Ке, д=1 она совпадает с поверхностью соответствующего тела). Выберем на каждой из вспомогательных поверхностей Бе, д конечную совокупность точек {МП д} д МП
и в каждой точке
п,д разместим пару независимых вспомогательных элементарных электрических диполей с моментами
n,q _ rin,^-^n,q
РтІ = РтІ
ТІ
n,q _ _ r.n,q-^n,q ет2 ,
Рт2 = = Рт2
ориентированными
выбранных
n,q ¿n,q
Єт2 :
области De, ке=ю(єеце) ; R
'MMnq ' n,q „n,q
пд х ЁЦд =-Пд X Ё'д, д = 1,2,..., 0, ] = 1,2,..., Ьд , (5)
где Пд, Ё3ед и Ё0д - значения векторов нормали и
электрических компонент рассеянного (4) и возбуждающего полей в точке ] на поверхности тела с номером д. Решение системы (5) определяется путем минимизации функционала
е I |2
Ф = ЁЁ\П]д х (Ё,д + Ё0,д)| . (6)
д=1 ] =1
Задача минимизации функционала (6) решается методом сопряженных градиентов.
После решения задачи минимизации (определения неизвестных дипольных моментов
рП1д, Ртп2д (д = Ш, п = )
необходимые характеристики рассеянного поля определяются из (4). В частности, для компонент рассеянного поля в дальней зоне имеем
Ёе,е(М) = ( /ее)1/2 Ие ф(М) =
= (ехр(ВД / ^Я) Ве (е, Ф) + 0(Я-2),
вдоль единичных направлений ет1 в плоскости, касательной к Бед в точке Мпд , и излучающих в однородную среду с параметрами ее, це. Представим неизвестное рассеянное поле {Ёе, не} в
Ве в виде суммы полей введенных вспомогательных диполей:
- т 2 Ыд
Ёе (М) = — ££Ух (УхП п,д),
^е д=1 п=1
1 е
Не (М) =— ХЕУхИ п,д ,
М е д=1 п=1 П п,д = Т е (М , Мп,д ) РГ ,
РП,д = рП1,деТП1,д + рП^де?Т2д, М е ве. (4)
Здесь Те (М,М„,д ) = ехр(1кеЯММп д ) / 4пЯММп,д - фундаментальное решение уравнения Гельмгольца для
£е,Ф (M) = -(( / Єе )І/2 яе,е (M ) = = (exp(/keR) / keR)Dy (Є, ф) + 0(R-2),
(7)
расстояние от точки
М„д до точки М в Ве, Рт1,д, Рт2д (д = 1, е, п = 1, N) -
неизвестные комплексные постоянные (дипольные моменты); - число точек размещения диполей на
вспомогательной поверхности Бед.
Поле (4) удовлетворяет уравнениям Максвелла (1) и условиям излучения (3) в области Ве. Для того, чтобы удовлетворить граничным условиям (2), необходимо соответствующим образом выбрать значения
дипольных моментов Ртп1д , Рп2 (д = 1, е, П = 1, Ыд ).
Используем для этого метод коллокаций. Пусть М]
( = \,2,...,Ьд) - точки коллокации на поверхности Бд, Ьд - число точек коллокации на Бд. Тогда для определения неизвестных Ртп1,д , Р-П2д (д = 1, е, П = 1, N) получим следующую систему линейных алгебраических уравнений с комплексной матрицей размерности:
( е А( е А
2Х Ьд х 2Х ^д V д=1 / V д=1
где компоненты диаграммы рассеяния Дз(9,ф) и -Оф(9,ф) определяются выражениями
Q Nq
D9 (9, ф) =XÍG»,? (9, Ф) Х
q=1 n=1
x{(cos9 cosфcosajn’q + cos9 sinфcosPjn’q -
- sin 9 cos Yjn’q)pTn1q +
+(cos9 cosфcos aj q + cos9 sinфcos pn’q -
- sin 9 cos y2q)Pa },
Q Nq
Ар (9, ф) = XÉGn,q (9, ф) Х
q=1 n=1
x{(- sin ф cos ajn’q + cos ф cos ) p”[q +
+(- sin ф cos a n,q + cos ф cos PJ q) p^ },
(8)
в которых cos a”,q ,cos Pj”,q ,cos yj”,q и cos a2’q ,cos el’q, cos y2’q - направляющие косинусы единичных векторов eT”j,q и eT”^q; xnq, yn,q, z„iq - декартовы координаты
точки Mnq q; 9 и ф - общепринятые угловые сферические координаты точки наблюдения M.
Контроль точности решения осуществляется путем вычисления относительного значения функционала (6) на сетке точек коллокации и точек, промежуточных к ним, выбираемых на поверхностях Sq всех тел, входящих в систему
Д = (Ф'/ Ф 0)І/2,
0 ¿q I . ^ .
Ф0 =SS Ч Х 4q\
q=i j=i
(9)
где Ф' - значение функционала (6) на указанной выше
259
2
совокупности точек, а Фо - значение соответствующей нормы падающего поля на этой же совокупности точек, L'q - число точек коллокации и промежуточных точек на поверхности рассеивателя с номером q.
Таким образом, построено некоторое решение граничной задачи (1) - (3). В соответствии с общей идеологией метода дискретных источников [4], для того чтобы решение (4) действительно являлось приближенным решением граничной задачи (1) - (3), достаточно доказать полноту системы функций
х 1 (M,Ыпл) = пх V х(VX(T(M,Ыпл)Cq))j ,
q=1, ... , Q, n=1, 2, ..., m=1, 2, (10)
в пространстве L2(Sj) x ...x L\(Sq ), где LL2( Sq) -
пространство векторных функций, лежащих в касательной к Sq плоскости, компоненты которых принадлежат L2(Sq). Основываясь на технике, развитой в [4] для одиночных идеально проводящих тел, нетрудно показать, что это имеет место, если вспомогательные поверхности Se q являются нерезонансными, а множе-
( ~iNq
ства точек {Mn q} ^ всюду плотными на соответствующих Se>q.
ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ
На основании изложенного выше метода создана программа для расчета компонент рассеянного поля и контроля точности полученного решения. Максимальное число тел, предусмотренное программой, равно пяти. Каждое из тел может быть либо трехосным эллипсоидом, либо конечным цилиндром с эллиптическим поперечным сечением (с плоскими или скругленными торцами), либо параллелепипедом; все острые кромки и вершины предполагаются скругленными. Сочетания выделенных геометрических форм в структуре могут быть произвольными
Помимо типов геометрии рассеивателей, входными величинами программы являются координаты точек Oq, определяющие взаимное расположение рассеивателей в структуре, величины характерных геометрических параметров (в длинах волн) выбранных
типов рассеивателей, возбуждающее поле {E), Hо },
параметры подобия Ke q, числа точек размещения диполей Nq и точек коллокации Lq для каждого из рассеивателей, образующих исследуемую структуру. Минимизация функционала (6) осуществляется методом сопряженных градиентов; итерационный процесс выполняется при условии, если относительное изменение функционала на каждой из десяти последних итераций не превышает о.оо1. При помощи данной программы выполнена серия вычислительных экспериментов, направленных как на выяснение влияния параметров метода на погрешность получаемого решения (на основе критерия невязки), так и на получение характеристик рассеяния конкретных структур. Некоторые результаты представлены ниже.
На рис. 2 приведены кривые, иллюстрирующие зависимость относительной нормы невязки граничных условий (9) от положения вспомогательных поверх-
ностей для структуры, состоящей из двух трехосных эллипсоидов с различными значениями полуосей. Центр первого эллипсоида совмещен с началом декартовой системы координат, центр второго эллипсоида размещен в точке с координатами (3.0628, 0, 0) т.е. на оси х. Полуоси эллипсоидов кга, кеЬ, кес ориентированы вдоль осей х, у, 2 соответственно и принимают значения кеа\=2.0, кеЬ 1=1.5, кес\=\.0 для первого эллипсоида и кеа2=1.0, кеЬ2=1.5, кес2=2.0 для второго. Легко видеть, что наименьшее расстояние между точками поверхностей взаимодействующих эллипсоидов равно 0.01Х, где X - длина волны возбуждающего поля во внешней среде. В качестве последнего выбрана линейно поляризованная плоская волна, распространяющаяся вдоль оси 2, электрический вектор Е0 которой ориентирован вдоль оси х.
Рис. 2. Зависимость нормы невязки граничных условий от положения вспомогательных поверхностей для структуры из двух трехосных эллипсоидов с полуосями ка = 2.0, кеЬ1 = 1.5, кес1 = 1.0 и кеа2 = 1.0, кеЬ2 = 1.5, кес2 = 2.0: 1 -Ке2 = 0.2; 2 - Ке2 = 0.6; 3 - Ке2 = 0.8
При решении задачи числа точек размещения диполей на вспомогательных поверхностях 5^,1 и 5е2 выбраны одинаковыми, М=Ж2=144, и распределены они в локальных сферических системах координат с центром в центрах соответствующих эллипсоидов следующим образом: в каждом из полусечений ф=соп81, общее число которых равно двенадцати, выбраны двенадцать точек размещения диполей равномерно по углу 9. Количества точек коллокации на поверхности каждого из эллипсоидов также выбраны одинаковыми: ¿1=£2=288. Число полусечений ф=сош1, в которых размещены точки коллокации, выбрано в два раза большим, чем при размещении диполей; расположение точек коллокации по углу 9 сохранено таким же, как для диполей.
Кривые 1, 2 и 3 на рис. 2 характеризуют зависимость нормы невязки граничных условий (9) от положения вспомогательной поверхности 5е,1 первого эллипсоида, характеризуемого значением параметра Кед, при различных фиксированных положениях вспомогательной поверхности 5е2 второго эллипсоида, характеризуемого значением параметра Ке2. Как видно из рисунка, при различных значениях Ке,2 наименьшее значение нормы невязки граничных условий
достигается при значениях параметра Кг\, лежащих в интервале 0.4 < Кг1 < 0.7. Аналогичные результаты получаются и для Ке,2, если в рассматриваемой задаче исследовать зависимость нормы невязки граничных условий от положения вспомогательной поверхности Бе2 при фиксированных положениях вспомогательной поверхности .
Подобные исследования проводились для структур, составленных из трех, четырех и пяти тел различной формы. Эти исследования показывают, что оптимальные значения параметров Ке,д (д = 1, Q), обеспечивающие минимальные значения нормы невязки граничных условий на совокупности идеально проводящих тел, лежат в интервале 0.4< Кед<0.7. Обычно без всякого дополнительного анализа достаточно выбрать Ке, д=0.5. Отметим, что установленный интервал оптимальных значений параметров подобия для совокупности идеально проводящих тел совпадает с интервалом оптимальных значений параметра подобия для одиночного рассеивателя [11].
Рис. 3 демонстрирует результаты сравнения сечений обратного рассеяния структуры в виде двух соприкасающихся сфер с соотношением радиусов 1.0:0.7, рассчитанных предлагаемым методом, с такими же результатами, приведенными в работе [12]. Структура размещена таким образом, что центры сфер расположены на оси г. Плоская волна падает на структуру вдоль оси г. По отношению к направлению прихода волны первой расположена сфера меньшего радиуса. По оси абсцисс на рис. 3 отложен волновой радиус большей сферы. При получении обсуждаемых результатов число точек коллокации было выбрано равным числу точек размещения диполей: Ь\=Ь2=Ы\=Ы2=128. Алгоритм размещения тех и других для обеих сфер был выбран одинаковым (в локальных системах координат): в каждом из шестнадцати полусечений ф=соп81, отстоящих друг от друга на угловое расстояние Дф=22.5°, равномерно по углу 9 располагалось восемь точек коллокации (точек размещения диполей).
Рис. 3. Сравнение сечений обратного рассеяния структуры в виде двух соприкасающихся сфер с соотношением радиусов 1.0:0.7, рассчитанных различными методами: 1 - результаты работы [12]; 2 - данный метод
Как показывают приведенные кривые, имеет место очень хорошее совпадение сравниваемых результатов. Небольшие наблюдаемые отличия можно объяс-
нить ошибками при графическом съеме информации с рисунка статьи [12], а также возрастанием ошибки наших расчетов (параметры метода фиксированы) при возрастании величины кеа; при кеа>3.0 норма невязки граничных условий (9) превышает 8%.
Рис. 4 и 5 характеризуют влияние электромагнитного взаимодействия между телами на бистатические сечения рассеяния для структур двух типов. Обе структуры составлены из двух трехосных эллипсоидов. Качественное отличие одной структуры от другой заключается в том, что первая структура состоит из двух одинаковых трехосных эллипсоидов, а вторая из двух различных трехосных эллипсоидов.
а0е/х2, дБ
0, град
Рис. 4. Бистатические сечения рассеяния структуры, составленной из двух одинаковых трехосных эллипсоидов с полуосями кеа1=кеа2=2.0, кеЬ1=кеЬ2=1.5, кес1=кес2=1.0: 1 - без и 2 -с учетом взаимодействия
О00/Х2,дБ
0, град
Рис. 5. Бистатические сечения рассеяния структуры, составленной из двух трехосных эллипсоидов с полуосями кеа1=2.0, кеЬ[=1.5, кес[=1.0 и кеа2=1.0, кеЬ2=1.5, кес2=2.0: 1 -без и 2 - с учетом взаимодействия
Для обеих структур полуоси эллипсоидов кеа, кеЬ, кес ориентированы вдоль осей х, у, 2 соответственно; структуры возбуждаются линейно поляризованной плоской волной, распространяющейся вдоль оси 2, с вектором Е0 , ориентированным вдоль оси х. В первой структуре эллипсоиды характеризуются полуосями кеа1= кеа2=2.0, кеЬ1= кеЬ2=1.5, кес1= кес2=1.0; центр одного из них размещен в точке (0, 0, 0), а другого в точке (4.0628, 0, 0) декартовой системы координат. Во второй структуре центр первого эллипсоида с полуосями кеа1=2.0, кеЬ1=1.5, кес1=1.0 также размещен в
начале декартовой системы координат, а центр второго эллипсоида с полуосями кеа2=1.0, кеЬ2=1.5, кее2=2.0 - в точке с координатами (3.0628, 0, 0). Таким образом, в обеих структурах центры эллипсоидов расположены на оси х, а наименьшее расстояние между точками поверхностей взаимодействующих эллипсоидов равно 0.01Х.
Параметры метода для каждого эллипсоида обеих структур выбраны одинаковыми: М=Ж2=168,
¿і=£2=336. В локальных системах координат с центрами в центрах эллипсоидов точки размещения диполей и точки коллокации распределены следующим образом. В каждом из четырнадцати полусечений Ф=сош1, отстоящих одно от другого на угловое расстояние Дф=25.7°, равномерно по углу 9 выбрано двенадцать точек размещения диполей. Для точек коллокации алгоритм их расположения по углу 9 выбран таким же, как для точек размещения диполей, но выбираются они как в полусечениях ф=сош1, определенных для точек размещения диполей, так и посередине межу ними. При данных параметрах метода для обеих структур значение невязки (9) не превышает 25%. Кривые 1 и 2 на рис. 4 и 5 - это бистатические сечения рассеяния описанных выше структур в Е-плоскости (плоскость хг). Кривые 1 - это результат
суперпозиции двух решений, одно из которых - решение задачи рассеяния на первом эллипсоиде структуры в отсутствие второго эллипсоида, а второе - решение задачи рассеяния на втором эллипсоиде в отсутствие первого. Очевидно, что такой подход не учитывает электромагнитного взаимодействия рассеивателей. Кривые 2 - результат решения задачи на двух телах предлагаемым методом в полной электродинамической постановке. Сравнение кривых 1 и 2 позволяет оценить влияние электромагнитного взаимодействия тел структуры на бистатические сечения рассеяния.
Интересно отметить (см. рис. 5), что наибольшие различия в диаграммах рассеяния проявились для структуры, состоящей из эллипсоидов с различными геометрическими параметрами, в областях, прилегающих к направлению бокового рассеяния (9«90°). Полученные результаты свидетельствуют о том, что, если рассеивающая структура содержит тела, расстояние между которыми много меньше длины волны возбуждающего поля, корректный анализ процессов рассеяния может быть выполнен только путем решения соответствующей граничной задачи в полной электродинамической постановке.
ЛИТЕРАТУРА
1. Cooray M.F.R., Ciric I.R. // Comput. Phys. Commun. 1991.V. 68. No. 1-3. P. 279.
2. YeeK., Chen J. // IEEE Trans. Antennas and Propag. 1997. V. AP-45. No.3. P. 354.
3. Hwu S. U., Wilton D.R., Rao S.M. // Procced of the Antennas and Propag. AP-S Int. Symp. Syracuse, N.Y., June 6-10, 1988. V.2. P. 890.
4. Еремин Ю. А., Свешников А. Г. Метод дискретных источников в задачах электромагнитной дифракции. М.: Изд-во Моск. ун-та, 1992.
5. Дмитренко А.Г., Мукомолов А.И. // РЭ. 1988. Т. 33. No.3. С. 449.
6. Дмитренко А.Г., Мукомолов А.И. // РЭ. 1990. Т. 35. No.2. С. 438.
7. Дмитренко А.Г., Мукомолов А.И. // РЭ. 1991. Т. 36. No.5. С. 1032.
8. Дмитренко А.Г., Мукомолов А.И. // РЭ. 1995. Т. 40. No.6. С. 875.
9. Дмитренко А.Г., Мукомолов А.И., Фисанов В.В. // РЭ. 1998. Т. 43. №8. С. 910.
10. Дмитренко А.Г., Корогодов С.В. // РЭ. 1998. Т. 43. №12. С. 1463.
11. Dmitrenko A. / Proceed of the 3rd Workshop on Electromag. and Light Scattering: Theory and Applications, Bremen, Germany, March 16 - 17, 1998. P. 59.
12. Свешников А.Г., Еремин Ю.А. // Численные методы электродинамики. М.: Изд-во Моск. ун-та, 1980. Вып. 4. С. 3.
Статья представлена кафедрой исследования операций факультета прикладной математики и кибернетики Томского государственного университета, поступила в научную редакцию 7 мая 2003 г.
