Фасеточная модель в задачах рассеяния электромагнитных волн на телах с импедансной поверхностью Текст научной статьи по специальности «Математика»

УДК 537.874.6

Н.Н. Олюнин

Московский физико-технический институт (государственный университет) Радиотехнический институт им. акад. А.Л. Минца РАН

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

поверхностью

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

Ключевые слова: рассеяние, фасеточная модель, импедансная поверхность.

I. Введение

II. Отражение от треугольного экрана

Для расчёта отражения электромагнитных волн от больших (в масштабе длины волны) тел можно использовать асимптотические методы теории дифракции: метод физической оптики [1] и метод краевых волн [2]. При этом, как правило, можно разделить поверхность тела на отдельные фрагменты, для которых в приближении этих методов получаются аналитические выражения для рассеянного поля. Удобной аппроксимацией поверхности тела является фасеточная модель. Согласно этой модели гладкая поверхность тела аппроксимируется треугольными площадками (фасетами), а изломы поверхности — прямыми рёбрами. Поле, рассеянное каждой фасетой и каждым ребром, рассчитывается отдельно, а затем эти поля суммируются с учётом взаимного расположения фасет и рёбер. Модель фасет требует значительных вычислительных ресурсов при расчёте отражения от протяженных тел сложной формы. Тем не менее эта модель часто используется на практике. Это связано прежде всего с возможностями синтеза фасеточного представления сложных объектов в современных системах автоматизированного проектирования [3].

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

:Гг:

Г =

Г =

где с — скорость света, Е и Н — векторы электромагнитного поля, п — нормаль к поверхности тела. Затем эти токи принимаются за источники рассеянного поля. Далее рассчитываются электрический и магнитный векторные потенциалы Ае и Ат путём интегрирования по поверхности тела:

Ае

1

- Ф

с

АкК

к

-Г-

ЧБ.

(1)

где к

волновое число

я

расстояние

от точки на поверхности тела до точки наблюдения. Вектор Е рассеянного поля выражается через векторные потенциалы следующим образом [4]:

Ё = — \^гасісііуАе + к2Ае) — тоіАт. (2) гк

Следует различать две поляризации падающей волны: Е-поляризация, когда вектор Е лежит в плоскости отражения, и Н-поляризация, когда в плоскости отражения лежит вектор Н.

Обозначим через р0 вектор поляризации падающей волны, через т0 и т — единичные векторы направлений падающей и отражённой волн, через р — единичный вектор, перпендикулярный т, через р = т — т0 — вектор рассеяния, а через КЕ и Ян — коэффициенты отражения от материала экрана для Е — и Н-поляризаций. Если поверхность идеально проводящая, то Яе = — 1, Ян = 1. В случае отражения плоской волны от поверхности с комплексной диэлектрической проницаемостью е коэффициенты отражения рассчитываются по формулам Френеля [4]. В случае отражения от поверхности с покрытием коэффициенты отражения рассчитываются по соотношениям, приведённым в [5]. Эти коэффициенты используются здесь, так как согласно приближению физической оптики предполагается, что каждая точка поверхности отражает падающую волну так, как её отражала бы касательная в данной точке плоскость. Использование приближения физической оптики в задачах рассеяния на телах с неидеальной проводимостью или покрытием правомерно только в случаях, когда толщина скин-слоя или покрытия мала по сравнению с размерами и радиусами кривизны тела.

Пусть экран образован двумя векторами и у2, имеющими начало в одной точке. Введём систему координат (X. У. Z), связанную с экраном (рис. 1). Ось У совпадает по направлению с нормалью к экрану п, векторы п и т0 лежат в плоскости ХУ, ось Z лежит в плоскости экрана и совпадает по направлению с векторным произведением [т0.п]. Пусть ех, еу и в2 — орты, соответствующие осям X, У и Z.

Рис. 1. Отражение от треугольного экрана

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

^ к РікЕ

где

1-Е

р = (роА )рЕ + ([Шо,ро],е2 )ри, (3)

1 + ЯЕ , ч 1 — ЯЕ ,

- {туЄ2-т2Єу)-\------— (?в0 ,п)е2

2

рп

1 — Ян 2

2

1 + Ян 2

ех +

(то,и)(тхЄу — Шувх),

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

N = фе-ікф,Ч8 = — аЬ

ае

іЬ

— Ье1

аЬ

1

где а = (—кр.щ), Ь = (—кр.щ), с = \ [щ р2]|. В случае идеальной проводимости выражение (3) приводится к следующему виду:

р = (т о .п)ро — (ро .п) т о.

III. Рассеяние на изломе поверхности

Для решения задачи рассеяния на изломе поверхности используется метод краевых волн, подробно описанный в [2]. Согласно этому методу, точное решение за-

дачи представляется в виде суммы равномерной и неравномерной частей. Равномерная часть соответствует решению, найденному в приближении физической оптики. Неравномерная часть является поправкой, обусловленной изломом поверхности. Асимптотический анализ показывает, что для расчёта поля в дальней зоне можно ввести неравномерные (рёберные) токи, являющиеся источниками неравномерной части поля.

Для задачи рассеяния на изломе поверхности эталонной является задача рассеяния плоской волны на бесконечном клине.

Введём цилиндрическую систему координат (р,р,г), связанную с клином (рис. 2). Ось г направим вдоль ребра, то есть перпендикулярно плоскости рис., плоскость р = 0 совместим с одной из граней клина, а плоскость р = а — с другой. Пусть на клин падает плоская волна

Е?

—ікр соб(^-^о)

или плоская волна

Н = е-ікр соъ(<р—<роо)

(4)

(5)

где ро — полярная координата направления на источник.

ном клине

Неравномерная часть поля, рассеянного клином, при кр ^ 0 имеет вид цилиндрических волн, расходящихся от ребра клина [2]:

ег(кр+п/4) Л(кр+п/4)

Ех = !1—==- и Н2 = д1

у/2тгкр

у/2тгкр

для волн (4) и (5) соответственно, где функции f1 и д1 зависят от углов р, р0 и а,

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

д = д — д

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

Рассмотрим теперь наклонное падение плоской волны на клин конечной длины, ребро которого совпадает с отрезком [—^/2; р/2] на оси г. Пусть р и р0 — полярные координаты направлений на источник и приёмник соответственно, а 7 — угол между направлением падающей волны и ребром (рис. 3). Введём прямоугольную систему координат (X, У, Z). Ось Z совпадает с осью г цилиндрической системы координат, ось X лежит в плоскости р = 0. Обозначим через р0 вектор поляризации падающей волны, через т0 и т — единичные векторы в направлениях падающей и отражённой волн, через рт — единичный вектор, перпендикулярный т, через р = т — т0 — вектор рассеяния.

Рис. 3. Наклонное падение плоской волны на клин

Будем предполагать, что на клине конечной длины наводятся такие же неравномерные токи, как и на бесконечном. Применяем формулу (1) для конечного и для бесконечного клина. Затем, используя формулу (2) и выражение для неравномерной части поля, рассеянного бесконечным клином при наклонном падении на него плоской волны, приведённое в [2], получим выражение для неравномерной части по-

0

ля, рассеянного конечным клином. Проекция напряжённости рассеянного поля на вектор р в дальней зоне равна

^ 1 е^кК

Б = (М&,р)—8шс(кдг(1/2)—-,

2п Я

где

M

f—g1 sin p0 sin p g1 cos p0 sin p —f1 cos pctgY\

V

g sin p0 cos p g1 cos p0 cos p —f1 sin pctgY n n fi

/

вектор рт0 задан координатами в системе (X, У, Z), а Я — расстояние от начала отсчёта системы (X, У, Z) до точки наблюдения.

IV. Выражения для функций / и 9

Для идеально проводящего клина выражения для функций / и д приведены в [2].

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

du . л ~—Ь гк sm в і 2и, dn

О.

(6)

^=0,а

гДе &1,2 — комплексные постоянные, ^ -дифференцирование в направлении нормали, внешней по отношению к клину, а и — либо Ех, либо Нх. Точное решение задачи рассеяния плоской волны на бесконечном клине при таких граничных условиях подробно изложено в [7] (константы в^2 соответствуют в1,2 в [7]). Согласно этому решению

д \ = (фо(Ф — Р — п)а(Ф — р — п) —

— Ф0(Ф — p + п)а(Ф — p + п))х

х(ф0(Ф — po))

i

(7)

где

Ф0(л) = Ф^Д )Ф(л — 2Ф,л2),

Ф^,6) = фф (Л+Ф+п/2—9)фф(л+Ф—п/2+9),

a(z)

COS[fl(Ф - Po)]

sin f!Z — sin[i^ — p0)]

1

2Ф ’

фф — функция Малюжинца. Для вычисления функции f в формулу (7) нужно подставлять константы в1}2, соответствующие падающей волне (4), а для вычисления д — соответствующие падающей волне (5).

Функция Малюжинца — это специальная функция, введённая Г.Д. Малюжин-цем. Ссылки на методы вычисления этой функции приведены в [6].

V. Выражения для функций /0 и д0

Для идеально проводящего клина выражения для функций ^ и д0 также приведены в [2].

Для клина с импедансными граничными условиями (6) функции f0 и д0 получены по формулам (1) и (2) путём интегрирования по полуплоскости. Для них получается следующее выражение:

f0

1 (І — Ri) sin p0 — (І + Ri) sin p

2 cos p0 + cos p

Si+

+ -((l-i?2) sill(Q'-p0)-(l + -R2) sill(Q'-p))x

x(cos(a — p0) + cos(a — p))-162, (8)

где $i,2 = 1, если падающая волна освещает первую (р = 0) или вторую (р = а) грань соответственно, и ^,2 = 0 в противном случае, R12 — коэффициенты отражения от первой и второй граней, соответствующие граничным условиям (6) и равные

cos Xl,2 — sin 01,2

cos Xi,2 + sin 6*1,2

где х1>2 — углы падения, то есть углы между внутренними нормалями к поверхностям и направлением падающей волны. Для вычисления функции f0 в формулу (8) нужно подставлять константы в1,2, соответствующие падающей волне (4), а для вычисления д0 — соответствующие падающей волне (5).

VI. Импедансные граничные условия

0

g

2

Импедансными граничными условиями являются граничные условия Леонто-вича [4], приближённо верные при больших абсолютных значениях диэлектрической проницаемости е. Граничным условиям Леонтовича соответствуют коэффициенты отражения

R

E

R

cos х + Vs' Н cos X' + l/v^

(x — угол падения). Таким образом, при использовании граничных условий Леон-товича комплексные константы из (6) будут следующими: для падающей волны (4) sin в = \/е, а для падающей волны (5) sin 9=1/ у/е.

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

cos \ - (iy/e)/(tg(y/ekcl))

cos х ~ l/v^

R

E =

Rh

«•"S \ + (iy/e)/(tg(y/£kd))'' COSX - (tg(y/£kd))/(iy/£)

cosx + (bg(V£kd))/(iy/e)'

Таким образом, комплексные константы из (6) в случае такой поверхности будут следующими: для падающей волны (4) sin 9 = iy/e/tg(y/ekd), а для падающей волны (5) sin в = tg(y/skd)/(iy/s).

VII. Сравнение результатов расчёта с эталонными значениями

Для сравнения результатов расчёта, полученных с помощью фасеточной модели, с эталонными значениями, необходимо ввести несколько обозначений. Введём прямоугольную систему координат (X, Y, Z) и сферическую систему координат (r, 0, р): зенитный угол 0 отсчитывается от оси Z, азимутальный угол р — от оси X. Пусть т — единичный вектор в направлении отражённой волны с угловыми координатами 0 и р. Углам 0 и р поставим в соответствие единичные векторы вв и ev, перпендикулярные т (рис. 4). Вектор вв лежит

в плоскости, образованной осью Z и вектором т, и направлен в сторону увеличения угла в. Вектор е^ равен [т,ве]. Разложим напряжённость поля по этим векторам: Е = Еовв + Е^е^. В дальней зоне для падающего поля Еі и рассеянного поля Е5 справедливо соотношение

Es

ES

S

щ,

El

kR

R

где S — двумерная матрица рассеяния,

13

Sii S12

S2i S22

(9)

Элементы матрицы рассеяния зависят от направлений падающей и отражённой волн.

Рис. 4. Обозначения для сравнения результатов расчёта с эталонными значениями

Будем полагать, что векторы направлений падающей и отражённой волн находятся в плоскости XZ. Тогда направление на приёмник определяется только одним углом в.

VIII. Сравнение с эталонными значениями для шара

Шар является одним из немногих геометрических объектов, для которых известно аналитическое решение задачи рассеяния.

Точное решение для шара рассчитывалось с использованием алгоритма, приведённого в [8], — путём суммирования рядов Ми, содержащих коэффициенты отражения для сферических гармоник.

На рис. 5 изображна зависимость поля, рассеянного шаром, от угла в, определяющего направление отражённой волны. Направление на источник совпадает с осью Z. Радиус шара а выбран так, что ка = 10. Диэлектрическая проницаемость материала шара е = 1.5 + 2,5г. Для сравнения приведена аналогичная зависимость для идеально проводящего шара такого же радиуса. Цифры, которыми отмечены кривые на

угол в

Рис. 5. Сравнение результатов расчёта с помощью фасеточной модели с точными значениями для шара

рис. 5, означают следующее: 1 — расчёт с использованием модели фасет, 2 — расчёт по формулам Ми для шара с конечной проводимостью, 3 — расчёт по формулам Ми для идеально проводящего шара.

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

IX. Сравнение с эталонными значениями для цилиндра

Эталонные значения поля, рассеянного цилиндром, рассчитывались с помощью программы, доступной в Интернете по адресу http://wave-scattering.com (автор программы — А. Мороз). Эта программа представляет собой модификацию другой программы, автор которой М.И. Мищенко. Изменения в программе Мищенко описаны в [9]. Программа Мищенко доступна в Интернете по адресу

http://www.giss.nasa.gov/~crmim. Пояснения, необходимые для работы с программой Мищенко, приведены в [І0]. Вычислительная процедура в этой программе основана на методе T-матриц. Обзор этой методики приведён в [ІІ].

Теперь перейдём к рассеянию на цилиндре. Пусть центр симметрии цилиндра совпадает с началом отсчёта системы (X, Y, Z). На рис. 6 и 7 приведены зависимости абсолютных значений двух компонент матрицы (9) от угла 6 при обратном рассеянии на наклонённом цилиндре. Приёмник совмещён с источником (в

этом случае Б12 = Б21). Выбраны следующие значения параметров цилиндра: радиус Я = ЗА (Л — длина волны), длина цилиндра Ь = 6А, е = 2,25 + 10г. Ориентация оси цилиндра определяется угловыми координатами в = п/3 и р = п/6. Цифры, которыми отмечены кривые на рис. 6 и 7, означают следующее: 1 — расчёт с использованием модели фасет для цилиндра с конечной проводимостью, 2 — расчёт с помощью метода Т-матриц для цилиндра

с конечной проводимостью, 3 — расчёт с использованием модели фасет для идеально проводящего цилиндра, 4 — расчёт с использованием модели фасет для цилиндра с конечной проводимостью без учёта рассеяния на кромках (кривые 3 и 4 приведены для сравнения).

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

Рис. 6. Рассеяние на круговом цилиндре. Зависисмость Зц(в)

X. Заключение

Основным достоинством описаной методики вычисления рассеянного поля является её универсальность в отношении формы тела. Главным отличием от аналогичных методов, рассмотренных в других работах, является использование обобщённого метода краевых волн.

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

Рис. 7. Рассеяние на круговом цилиндре. Зависисмость Si2(d)

Литература

1. Потехин А.И. Некоторые задачи дифракции электромагнитных волн. — М.: Советское радио, 1948.

2. Уфимцев П.Я. Метод краевых волн в физической теории дифракции. — М.: Советское радио, 1962.

3. Борзов А.Б., Засовин Э.А., Соколов А.В., Сучков В.Б. Методы синтеза геометрических моделей сложных радиолокационных объектов // Электромагнитные волны и электронные системы. — 2004. — Т. 8, № 5.

4. Вайнштейн Л.А. Электромагнитные волны. — М.: Советское радио, 1957.

5. Бреховских Л.М. Волны в слоистых средах. — М.: Наука, 1973.

6. Osipov A. A simple approximation of the Maliuzhinets function for describing wedge diffraction / / IEEE Trans. on Antennas and Propagation. — 2005. — V. 53, N. 8. — P. 2773-2776.

7. Бабич В.М., Лялинов М.А., Грику-ров В.Э. Метод Зоммерфельда-Малюжин-ца в теории дифракции. — СПб.: СПБГУ, 2003.

8. Борен К., Хафмен Д. Поглощение и рассеяние света малыми частицами. — М.: Мир, 1986.

9. Moroz A. Improvement of Mishchenko’s T-matrix code for absorbing particles // Applied Optics. — 2005. -V. 44, N. 17. — P. 3604-3609.

10. Mishchenko M.I. Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation // Applied Optics. — 2000. — V. 39, N. 6. — P. 1026-1031.

11. Mishchenko M.I., Travis L.D., Mackowski D.W. T-Matrix computations of light scattering by nonspherical particles: a review // JQSRT. — 1996. — V. 55, N. 5. — P. 535-575.

Поступила в редакцию 13.01.2008.