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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ |

УДК 535.36

В. И. А л е х н о в и ч, К. И. Зайцев, В. Е. К а р а с и к, И. Н. Фокина

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА РАССЕЯНИЯ ЭЛЕКТРОМАГНИТНОГО ИЗЛУЧЕНИЯ НА ПРОВОДЯЩИХ ДИЭЛЕКТРИЧЕСКИХ ТЕЛАХ СЛОЖНОЙ ФОРМЫ

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

E-mail: [email protected]

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

разностей, дифракционный интеграл Релея-Зоммерфельда.

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

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

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

10 ... 20 длин волн от центра рассеивателя) путем решения уравнений Максвелла методом конечных разностей (finite-difference time-domain (FDTD)). Второй этап моделирования предполагает определения поля в дальней зоне путем вычисления дифракционного интеграла от распределения комплексной амплитуды поля на контуре, охватывающем расчетную область. Распределение поля в дальней зоне позволит определить дифференциальное сечение рассеяния, а также индикатрису рассеяния.

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

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

Запишем уравнения Максвелла в разностной форме

где Е — вектор напряженности электрического поля; Н — вектор напряженности магнитного поля; е — относительная диэлектрическая проницаемость среды; ц — относительная магнитная проницаемость

Лоигсе — вектор плотности электрического тока; а* — эквивалентные магнитные потери; а — электрическая проводимость среды [2].

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

(1)

среды; Ms

source

вектор эквивалентной плотности магнитного тока;

(2)

Так как структура и свойства объекта неизменны в направлении оси OZ, а рассматриваемые среды являются линейными, то рассеянное на объекте излучение будет иметь ТМ поляризацию:

'Ех = 0\ (Их = 0

И = ( И = 01. (3)

\И2 = 0

E =

Ey = 0

Ez = 0,

Принятые допущения позволяют записать уравнения Максвелла в следующем виде:

дИх dt

dHy

1 V

dEz

1 V

_Z = 1 dt е

dt DE.

dy dEz

— (Msourcex + a*Hx)

dx

- M

+ o*Hy

(4)

дНУ dHx _ ( т + E ч

(Jsourcez + " Ez )

dx

dy

Падающая i

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

продолжительность процесса численного моделирования. На рис. 1 приведена схема моделируемого пространства при принятых допущениях.

Расчетная область дискретизируется по пространственным и временной координатам в соответствии с алгоритмом Yee [3]. Шаг пространственной дискретизации Ах = Ду для обеспечения корректности моделирования выбирается из условия:

Amin

Ах = Ау <

20

(5)

Рис. 1. Схема моделируемого пространства

где Amin = c/(fmax^emax) — минимальная длина волны излучения; fmax — максимальная временная частота излучения; emax — максимальная диэлектрическая проницаемость в моделируемом пространстве.

Дискретизация по времени определяется в соответствии с условием устойчивости

y

Куранта-Фридрихса-Леви [2]:

At <

1

Ax

+

1

АУ

Ax

cT2'

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

Для введения поля в расчетную область, а также для нормирования параметров рассеянной электромагнитной волны воспользуемся методом полного и рассеянного полей (total field / scatter field (TF/SF)) [2]. В соответствии с указанным методом в расчетной области выделяется замкнутый контур, охватывающий рассеиватель (рис. 2).

Внутри TF/SF контура существует как падающее на частицу поле, так и поле, рассеянное частицей:

Еtotal — Einc + Е scatj

—total — —inc + Н;

(7)

sca^

где^щеи Hiric — падающее поле, £scat и Hiscat — рассеянное поле, Etotai и Htotal — полное поле. Вне TF/SF контура присутствует только рассеянное поле.

Расчет полей осуществляется в двух модельных областях пространства одинаковых размеров, в одной из которых рассеиватель отсутствует, а в другой — присутствует. В пустой расчетной области задается источник электромагнитного поля, причем плоская волна вводится с помощью метода мягкого источника [4]. Электромагнитная волна из вспомогательной расчетной области вводится через TF/SF границу в пространство, содержащее рассеиватель. Рассчитывается полное поле после дифракции волны на частице. Параллельно с вычислением

Рис. 2. Иллюстрация к методу полного и рассеянного поля [2]

2

2

полного поля моделируется распространение плоской волны в пустой расчетной области. На границе ТБ/ЗБ контура осуществляется вычитание падающего поля (поля пустой расчетной области) из полного поля (поля расчетной области, содержащей рассеиватель), в результате чего вне ТБ/ЗБ контура мы имеем только рассеянное частицей поле.

В результате численного моделирования взаимодействия электромагнитной волны с объектом определяется зависимость векторов напряженностей электрического Е(£) и магнитного Н(£) полей для каждой точки расчетной области. От временных зависимостей полей можно перейти к комплексной амплитуде электрического ЕЕ (и)

и магнитного НН(и) полей, используя алгоритм быстрого фурье-преобразования.

Преобразование поля ближней зоны в поле дальней зоны. Окружим область ТБ/ЗБ контуром произвольной формы, и восстановим комплексные амплитуды векторов напряженностей электрического и магнитного полей на данном контуре, полагая, что излучение имеет частоту и = и':

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

Для вычисления поля в дальней зоне воспользуемся дифракционным интегралом Релея — Зоммерфельда, который применительно к излучению с ТМ поляризацией имеет вид:

здесь Са — контур, г — радиус вектор, соединяющий начало координат и точку анализа поля в дальней зоне, г' — радиус вектор, соединяющий начало координат и точку на контуре, пСа — единичный вектор-нормаль контура.

В (9) используется функция Грина цилиндрической волны, имеющая вид [2]:

(8)

(10)

(2)

где H0 ' (x) — функция Ханкеля 2-го рода нулевого порядка.

При k\r — r '\ —у то, можно переписать функцию Грина (10) в виде

=-jk| r-r

G (r, r ') =

. 3

j2 e

V8nk\r _ r'

i •

_ r ' \ 2

(11)

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

Ez (r) =

,-jkr gjn/4

fr f8nk

шрог' na x H (r') +kr' x na x E (r')

x e+jkr^ dC =

+

x

e-jkr ej^/4 Г \ шроУ • J (f) —

f8nk J

Ca

—kr' x M (r') • r

e+jkr^ dC, (12)

где р0 — магнитная постоянная; г' — орт-вектор оси 02; г — орт-

вектор, направленный на точку анализа; ,1 (г') и М (г') — эквивалентные электрический и магнитный токи.

Переходя в цилиндрическую систему координат, можно получить

зависимость интенсивности электрического поля на цилиндре радиу-

2

. Выражение для определения дифференциального

сом r = р: Ez (ф) сечения рассеяния частицы имеет вид[5]:

о (ф) = lim 2пр

р^-ж

Ez (ф)

Efc

(13)

где Ё™ — комплексная амплитуда падающего поля; р — расстояние от начала координат до точки анализа в дальней зоне.

Частица характеризуется полным сечением рассеяния, которое может быть найдено интегрированием дифференциального сечения в пределах угла 2п [1]:

-2п

°р =

o(y>)dy>.

(14)

Индикатриса рассеяния частицы связана с сечением рассеяния следующим образом [1]:

X М = (15)

Ор

Верификация алгоритма численного моделирования рассеяния путем решения тестовой задачи с известным аналитическим решением. Для верификации алгоритма была решена задача с из-

2

o

вестным аналитическим решением. Были найдены дифференциальное сечение рассеяния и индикатриса рассеяния частицы цилиндрической формы, которая имеет круглое сечение Z = Z'. Рассматривалась частица с диаметром О = 1,6А и диэлектрической проницаемостью е = 4. Поглощение частицы при решении тестовой задачи принимается равным нулю, поскольку рассматриваемый нами аналитический метод определения индикатрисы используется только для непроводящих сред.

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

Дифференциальное сечение рассеяния цилиндрической частицы, полученное путем аналитического расчета, приведено на рис. 3, б. Расчет осуществлялся в соответствии с известными выражениями (16)-(18) [1].

4 то 2

кпапсоэ (пц>)

п=0

'1, п = 0; 2, п = 0,

(к4а) Зп (коа) - щЗ'п (коа) Зп (к¿а)

° = k

ko

k —

""П

(16)

(17)

(18)

щЗп (каа) И(2) (коа) - пЗ (к^а) И(2) (коа) где а — радиус цилиндра; Зп (ж) — функция Бесселя п-го порядка; и(2) (ж) — функция Ханкеля 2-го рода п-го порядка; к^ = ш^ёфО — волновое число при распространении волны в рассеивателе, щ =

Г^оо „ опп

= а--характеристический импеданс материала частицы, по = 377 —

V е^

характеристический импеданс свободного пространства.

Рис. 3. Дифференциальное сечение рассеяния цилиндрической частицы с D = 1,6Л, е = 4

Угловая координата, градус

Рис. 4. Индикатриса рассеяния цилиндрической частицы

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

На рис. 4 приведена индикатриса рассеяния цилиндрической частицы, полученная в результате численного моделирования в соответствии с [5].

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

СПИСОК ЛИТЕРАТУРЫ

1. Исимару А. Распространение и рассеяние волн в случайно-неоднородных

средах. Т.Н. - М.: Мир, 1981. - 280 с.

2. H a o Y., M i 11 г a R. FDTD Modeling of Metamaterials. Theory and Applications.

Boston: Arteck House, 2009. - 379 p.

3. Taflove A., Hagness S. C. Computational Electrodynamics: The Finite-Difference Time-Domain Method. Boston: Arteck House, 2000.

4. Schneider J. B./ Understanding the Finite-Difference Time-Domain Method. Online book. 2012 URL http://www.eecs.wsu.edu/schneidj/ufdtd (последнее обращение 23.07.2012).

5. UmashankarK., Taflove A. A Novel Method to Analyse Electromagnetic Scatterinf of Complex Objects // IEEE Transactions on Electromagnetic Compatibility. - 1982. - V. EMC-24, no. 4. - P. 397-405.

Статья поступила в редакцию 27.07.2012