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

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ЭЛЕКТРОМАГНИТНЫХ ВОЛН В СРЕДЕ С ИСПОЛЬЗОВАНИЕМ МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ НА СМЕЩЕННЫХ СЕТКАХ Скубачевский А. А.

Московский физико-технический институт, http://mipt.ru 141700 Долгопрудный, Московская область, Российская Федерация

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

Представлена чл.ксрр. РАН и РАЕН И.Б. Петровым

Метод численного решения системы уравнений Максвелла конечными разностями на смещенных сетках — Finite Difference Time Domain (FDTD) второго порядка аппроксимации, известный также как Yee algorithm, был использован для создания пакета программного обеспечения для моделирования распространения электромагнитных волн в неграфической дисперсионной среде с условием устойчивости Куранта. Описаны основные преимущества этого метода. Получены разностные выражения для уравнений Максвелла для 2D и 3D случаев. Рассмотрены различные граничные условия — периодические, отражающие и поглощающие. Приведены численные уравнения для поглощающих граничных условий (PML). Неоднородная среда также входит в модель, показаны примеры работы созданной программы в ней. В качестве источников рассмотрены гармоничный источник, гауссиан и плоская волна. Плоская волна промоделирована с помощью техники Total Field-Scattered Field (TFSF). Тестовые расчеты сделаны для гармонических источников, для рупорной антенны и для плоской волны, распространяющейся через проводники различных проводимостей и через диэлектрики. Программа показала достаточную надежность в работе и широкие возможности для дальнейшего развития.

Ключевые слова: уравнения Максвелла, алгоритм Йи, электромагнитные волны, TM-, ТЕ-волны, алгоритм хорошо сопряженных слоев (PML), техника полного поля/рассеянного поля (TFSF) УДК 517.9

Содержание

1. введение (73)

2. модификация йи алгоритма (74)

2.1. Конечно-разностные аналоги уравнений Максвелла (74)

2.2. Граничные условия (76)

2.3. Источники излучения (77)

2.4. Результаты расчетов (78)

3. заключение (79) литература (79)

1. ВВЕДЕНИЕ

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

ряд численных методов решения уравнений Максвелла. Одним из самых распространенных таких методов является метод конечных разностей на смещенных сетках — Finite Difference Time Domain (FDTM) [1], также известный как Yee algorithm [2]. Преимуществами данного алгоритма является, кроме прочего, простая декомпозиция по данным (данные хранятся в естественном, понятном для каждого виде, с ними легко работать), наличие данных о значениях электрического и магнитного полей на каждом шаге по времени и в каждой точке пространства, а также простота построения радиолокационных портретов объектов. Несмотря на широкую применимость данного метода, до сих пор остается нерешенной задача создания солвера FDTD на произвольных сетках. В

СКУБАЧЕВСКИИ А.А.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

настоящей работе предлагается алгоритм численного получения нестационарных решений уравнений Максвелла для задачи распространения электромагнитных

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

2. МОДИФИКАЦИЯ ЙИ АЛГОРИТМА 2.1. Конечно-разностные аналоги уравнений Максвелла

Система уравнений Максвелла:

^ = _у х Е _ М,

дг

dD

= VxH - J,

И

И

дЕ 1 е Й 1 Z? — = —VxH--аЕ.

dt s s

---- Ну

У Ех ■©-

0,1,к) н„

дг

УБ = 0, УВ = 0.

где Е — напряженность электрического поля, D — электрическая индукция, Н

— напряженность магнитного поля, В

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

Б = ее0 Е,

В = мМоН,

где £ — диэлектрическая проницаемость, £0 — диэлектрическая проницаемость вакуума, ц

— магнитная проницаемость, ц0 — магнитная проницаемость вакуума.

1 = стЕ,

М = ст* Н,

где о — электрическая проводимость, о* — аналогичная магнитная проводимость. Таким образом,

дН 1 ъ ? 1 * Й

-=--УхЕ--ст Н,

дг

Рис. 1. Ячейка сетки FDTD алгоритма. В методе FDTD (Finite-Difference TimeDomain) электрическое и магнитное поля смещены на сетке друг относительно друга так, что каждая компонента одного поля окружена компонентами другого поля (см. рис. 1).

В простейшем одномерном случае на первом шаге по времени задается поле E от некоторого типа источника (к примеру, гармонического, гауссиана или плоской волны). На следующем шаге по времени рассчитывается поле H во всех ячейках сетки, причем каждая компонента поля H окружена двумя (т.к. одномерный случай) компонентами поля E. На третьем шаге по времени мы уже рассчитываем поле E, зная поле H на предыдущем шаге, и т.д. (см. рис. 2).

t=2it

Г=1.5Д1

t=At

t=0.5At

t=0

х=0 х=Дх

Рис. 2. Схема

х=2Дх х=ЗДх

интегрирования по времени.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Разностные уравнения [2] получим исходя из дискретизации пространственных и временных производных их конечно-разностными аналогами:

— (/8х, ]8у, к8г, п81) = -дх

8х

- + 0[(Ах)2

Е_

Е 1п+1/2 + е |п-1/2

'"х I/,1,к х I/, 1 ,к

х /,1,к

полей:

Е |п+1/2 = с | е |п-1/2 +

х I/, 1+1/2,к+1/2 ^а I/, 1+1/2,к+1/2 х I/, 1+1/2,к+1/2 ^

+СЬ /, 1 +1/2,к+1/2 П1+1,к+1/2 — Нг П1 к+1/2 —

—Ну |п, 1+1/2,к+1 + Н у П1+1/2,к

Б

а I/, 1,к

Б

4 I/, 1 ,к~

д цП+1/2 - «п-1/2 — (/8х, ]8 у, к8г, п8г) = -1-^ + 0[(А/)2 ].

8х

При такой дискретизации поле в правой части полевых уравнений будет считаться в момент времени п (производные по координате), а в левой — в моменты п-1/2 и п+1/2. Представим поле в точке п в следующем виде:

Важными для рассмотрения являются случаи ТМ и ТЕ-волны. Для них численные уравнения будут иметь следующий вид: ТМ-волна:

В конечном счете, имеем для компонент

Е |п+1/2 = с | е |п-1/2 +

I/-1/2,1+1/2 а 1/-1/2,1+1/2 ^ г I/-1/2,1+1/2 т +С4 |/-1/2,1+1/2 (Ну /, 1+1/2 — Ну |/-1,1+1/2 п +

+Н |/-1/2,1 П - Н | п-1/2,1+1),

Н х 1-1/2,1+1 Ба |/—1/2,1+1 Нх |П -1/2,1+1 +

+Б

4 I/-1/2,1+1

(Е Г

1/2,1+1/2

- Е.

Ну |п ,++1/2 Ба |/, 1+1/2 Нх |п 1+1/2

г I/-1/2,1+3/2 +

п +1/2),

+Д,

( Ег

- Е,

4 I/, 1I/+1/2,1+1/2 ^ I/-1/2,1+1/2

ТЕ-волна:

Е |п+1/2 = с | е |п-1/2 +

х I/, 1+1/2 а I/, 1+1/2 ^х I/, 1+1/2 т

+С4 |/, 1+1/2 (Н Г, 1+1 -Н2 |/, 1 п).

п +1/2).

Е

п+1/2

,= с.

Е

|п-1/2

у I/-1/2,1+1,к+1/2 _ а I/-1/2,1+1,к+1/2 ^у I/-1/2,1+1,к+1/2 +С4 |/-1,1+1,к+1/2 (Нх |п -1/2,1+1,к+1 -Нх |п -1/2,1+1,к -

-|п 1+1,к+1/2 + 1-1,1+1,к+1/2 X

|п+1/2 "г |/-1/2,1+1/2,к+1 "

| п-1/2

-а |/-1/2,1+1/2,к+1 Ег |/-1/2,1+1/2,к+1

+С4 |/-1/2,1+1/2,к+1 (Ну |1,1+1/2,к+1 Ну 1-1,1+1/2,к+1 -Нх / -1/2,1+1,к+1 + Нх |/ -1/2,1,к+1)'

Нх 1-1/2,1+1,к+1 Ба /-1/2,1+1,к+1 Нх 1-1/2,1+1,к+1 +

+Б

(Е |п+1/2 \i-1l

- Е.

п+1/2

4 I/-1/2,1+1,к+М у I/-1/2,1+1,к+3/2 1/-1/2,1+1,к+1/2

' |п+1/2 + е |п+1/2 )

г |/-1/2,1+3/2,к+1 + Ег |/-1/2,1+1/2,к+1 )'

Ну п,1+1/2, к +1 Ба |/, 1+1/2,к +1 Н у |п1+1/2,к +1

+Б

4 /, 1 +1/2,к+1

г+1/2

I/, 1 +1/2, к +3/2

(Е Г,

1/2

+ 1/2

+Е„

1/2,1 -1/2,к+1 I/-1/2,1+1/2,к +1

|„+1/2 . I/, 1 +1/2, к +1/2/'

где (размер ячейки сетки А):

с

а /,1,к

с

4 I/, 1,к

е. ., А

V /,1,к

|п+1/2 "у |/-1/2,1+1/2"

| п-1/2

"а |/-1/2,1+1/2 Еу |/-1/2,1+1/2

+с,

4 I/-1/2,1+1/2

(Н г ,+1 -Н

г I/, 1+1

п),

Н г / ,1+1/2 Ба |/, 1+1/2 Нг |п 1+1/2 +

+А (Е |п-Е, к,,, п +1/2 +

4 I/, 1+1/2

1+3/2

"г I/, 1+1/2

п+1/2).

+Е |п+1/2 -Е |

I/-1/2,1+1 \/+1/2,1+1

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

СКУБАЧЕВСКИИ А.А.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Условие сходимости схемы (условие Куранта):

111,

ст —- + — + —- < 1, VA? Al A?

1

где

с =

скорость света в среде.

Закон численной дисперсии:

1 . (а>Ы

-sin I-

с At I l

(

1 • fkyAy sin

Ay

//

1 . f к Ax —sin I —— Ax I l

,1 • f К Az + I — sin I —— Az I l

2.2. Граничные условия

В разработанной программе реализовано 3 типа граничных условий: отражающие, периодические и неотражающие.

Самыми простыми в реализации являются отражающие граничные условия. Самый их естественный вид — "идеальный проводник" (perfect electric conductor). Для их реализации достаточно задать на границе значение поля, равное нулю. Другой вариант отражающих граничных условий — имитация столкновения с такой же волной, но идущей в противоположном направлении: задание в граничных ячейках значений поля, равных значению в ячейках у границы. Также можно просто умножить значение поля в необходимых ячейках на очень малую величину (например, 10-10), но это, хоть и просто, но не физично, и может быть использовано только в тестовой программе. Последний способ моделирования отражения — задание высокой проводимости (см. о в уравнениях Максвелла) в необходимых ячейках.

Периодические граничные

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

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

Последний тип граничных условий — неотражающие. В работе используются неотражающие граничные условия, построенные на алгоритме Perfectly Matched Layer (PML) [2]. Основной идеей данных граничных условий является определенное задание проводимости на краях сетки, в слое толщины более одной ячейки (см. рис. 3). Преимуществом таких поглощающих граничных условий является то, что они верно работают при волнах любой поляризации, частоты и излучаемых любым источником (до 1994 года были также Absorbing Boundary conditions, но они работали только для плоских волн, и потому имели крайне ограниченное применение). Данный алгоритм основан на разделении каждой компоненты поля на две ортогональных друг другу. То есть вместо

РМЦСТ^.О;,.^,,^;)

h

рмм

PMLto^.ir^.O.O)

PML(oi,,(TIj.O.Q)

PMUff^.ir'i.e^.ll',) PMLlO.O.CT^.oj,) PMLICT^.trlj.ii^.o',)

Рис. 3. Граничные условия поглощения PML.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

системы из 6 уравнений, записанных ранее для 3-мерного случая имеем систему из 12 уравнений.

При реализации данного типа граничных условий есть нюансы, которые необходимо учесть, иначе они будут неверно работать (будет небольшое отражение): проводимость на краях области должна быть не постоянной во всех ячейках слоя PML, а возрастать от центра к границе сетки (например, как х3). Кроме того, минимизация величины отраженной волны (она может составлять как 10-2, так и 10-5 от величины падающей волны) требует отдельного исследования и тщательного подбора параметров.

В данном случае (как уже было сказано, каждая компонента поля разделяется на две) имеем систему уравнений в 3D:

5 a+с )Ey =

5 л ~+ С dt , E = ~ (^ + ),

5 dt % )E* = = | (Hxy + Hxz ),

5 N — + С dt x, E = (НУ + Hzz ),

5 л — + с dt x, E = 1 H + Hy),

5 E = = —дду(Hxy + Hz),

д * i — + су dt y J Hy

д „' i--+ с dt , J Hz - = д^ (+ E, ),

д „' l~dt + С J Hz = —д^ (Ey + Ez),

д „' i--+ с дt J Hx = дх (^ + ^ ),

д !-+ С* дг Х J Hzx-- = —дХ (E, + E, ),

Для случая TM и TE-волн можно довольно просто получить системы уравнений как частные случаи (см. [1]).

2.3.Источники излучения

В качестве источников в программе могут быть заданы гармонический источник, гауссиан и плоская волна. Задание первых двух не представляет никакой сложности (просто задание поля в определенных ячейках сетки). При моделировании плоской волны все обстоит немного сложнее. Ее можно промоделировать с помощью техники Total Field-Scattered Field (TFSF, [4, 5]; см. рис. 4). Основная сложность — необходимость задания источника, излучающего в одном направлении, а не симметрично в разных (в отличие от любого простого источника, заданного значением поля в определенных ячейках).

Это делается следующим образом. Расчетная область делится на области "Total field" - в которой есть как отраженное, так и поле от источника и "Scattered field" - поле только от источника. Граница между этими двумя областями и является "источником" плоской волны. Для реализации данной

Рис. 4. Области полного и отраженного полей.

78

70 СКУБАЧЕВСКИЙ А.А.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Рис. 5. Моделирование работы рупорной антенны. идеи задается поле во всех точках границы (например, как гауссиан), а в граничных точках области отраженного поля просто-напросто это поле вычитается из имеющихся значений. Поле Н, вычитаемое в области отраженного поля, если поле Е задано как гауссиан £(1), с учетом сдвига по времени и координате между полями:

- I'+?+т

2.4. Результаты расчетов

В данном разделе будут показаны результаты некоторых расчетов с помощью созданной программы.

В двумерном случае использовались следующие параметры расчетов: частота 5 ГГц; размер расчетной области 5х5 м (100x100 ячеек сетки; 5 см размер ячейки). На рис. 5 представлены результаты моделирования работы рупорной антенны. На рис. 6 можно наблюдать распространение излучения от двух источников, излучающих по гармоническому закону.

Далее на сетке 200х200 ячеек с размером ячейки 5 мм представлены примеры моделирования плоской волны (ТРЗР), проходящей через проводник (рис. 7),

®Н |Ц®:

'Д

|

-■йй»

Рис. 6. Распространение излучения от двух источников, излучающих по гармоническому закону.

Рис. 7. Моделирования плоской волны, проходящей через проводник.

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Рис. 8. Моделирования плоской волны, проходящей через диэлектрик.

Рис. 9. Падение плоской волны на проводник размером

на рис. 8 — через диэлектрик. на рис. 9а — через такой же проводник, но размера 1, 9Ь — проводимость меняется на 2 СМ/м, 9с -10 СМ/м, 9d - 100 СМ/м (волна падает на них так же, как и предыдущих случаях, поэтому на рисунках показано только то, что произошло с волной после прохождения проводника).

Амплитуда плоской волны 1 В/м; вид импульса можно видеть в одномерном разрезе на рисунке.

3. ЗАКЛЮЧЕНИЕ

Таким образом, в данном программном комплексе промоделировано

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

в 1 ячейку и пр

0.05, 2, 10, 100.

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

ЛИТЕРАТУРА

1. Taflove A, Hagness S. Computational electrodynamics: the finite-difference time-domain method (3d ed.) London, Artech House, 2005.

2. Yee Kane. Numerical solution of initial boundary value problems involving

s0

INFORMATION TECHNOLOGIES

Maxwell's equations in isotropic media. IEEE Transactions on Antennas and Propagation, 1966, 14 (3):302—307.

3. Rumpf RC, Garcia CR, Berry EA, Barton JH. Finite-Difference Frequency-Domain Algorithm for Modeling Electromagnetic Scattering from General Anisotropic Objects. PIERS B, 2014, 61:55-67.

4. Holland R, Williams J. Total-field versus scattered-field finite-difference Codes: A Comparative Assessment. IEEE Trans. Nuclear Science, 1983, 30(6):4583-4588, DOI: 10.1109/TNS.1983.4333175.

5. Berenger J-P. A perfectly matched layer for the absorption of electromagnetic waves. J. Computational Physus, 1994, 114(2):185-200.

Скубачевский Антон Александрович

студент

Московский физико-технический иснтитут, факультет управления и прикладной математики, кафедра информатики 141701 Долгопрудный, Моск. обл., Россия [email protected]

NUMERICAL SOLUTION OF MAXWELL EQUATIONS FOR MODELING OF ELECTROMAGNETIC WAVES PROPAGATION Anton A. Skubachevsky

Moscow Institute of Physics and Technology, http://mipt.ru

Faculty of Management and Applied Mathematics, Department of Computer Science 9, Institutskii per., 141701 Dolgoprudny, Moscow Region, Russian Federation [email protected]

Abstract. The method of numerical solution of Maxwell's equations FDTD second-order approximation. The brief review made in this area and possible applications. The advantages of this method. The derivation of difference equations used in this method. When covering a variety of boundary conditions and the characteristics of their modeling. For non-reflective boundary conditions presented differential equations used for their modeling. Examples of simulation of electromagnetic waves with different characteristics of the medium (simulated waves passing through the dielectric and conductor) and sources of radiation, including the horn antenna, the two-dimensional and three-dimensional cases. For clarity, also shows the propagation of waves in one-dimensional perspective. The radiation sources used Gaussian, harmonic source and a plane wave.

Keywords: Maxwell equations, Yee algorithm, electromagnetic waves, TM-wave, PML, TFSF. PACS: 02.60.Cb, 02.70.-c

Bibliography — 5 references Received 28.03.2016 RENSIT, 2016, 8(1):73-80_DOI: 10.17725/rensit.2016.08.073