УДК 517.9
А. A. Скубачевский, Н. И. Хохлов
Московский физико-технический институт (государственный университет)
Численное решение уравнений Максвелла для моделирования распространения электромагнитных
волн
В статье рассмотрен метод численного решения системы уравнений Максвелла FDTD второго порядка аппроксимации. Проведен краткий обзор сделанного в этой области и возможных приложений. Показаны преимущества данного метода. Представлен вывод разностных уравнений, используемых в методе. Освещены различные граничные условия и особенности их моделирования. Для неотражающих граничных условий представлены разностные уравнения, используемые для их моделирования. Приведены примеры моделирования электромагнитных волн при различных характеристиках среды (промоделировано прохождение волны через проводник и диэлектрик) и источниках излучения, в том числе рупорной антенны, в двумерном и трехмерном случаях. Для наглядности также показано распространение волн в одномерном разрезе. В качестве источников излучения использовались гауссиан, гармонический источник и плоская волна.
Ключевые слова: уравнения Максвелла, алгоритм Йи, электромагнитные волны, TM-волна, PML, TFSF.
A. A. Skubachevskii, N. I. Khokhlov
Moscow Institute of Physics and Technology (State University)
Simulation of electromagnetic waves propagation using the
FDTD method
The second order of approximation numerical method for solution of Maxwell equations in nondispersive medium called FDTD (finite difference time domain) also known as the Yee algorithm, is used to create the software package for modelling of electromagnetic waves propagation in medium (with the Courant stability condition). The total field/scattered field (TFSF) source is considered in detail. Inhomogeneous medium is also included in the model and there are some examples which show how the program works with it. The various boundary conditions (periodic, reflecting and absorbing) are considered. Perfectly matched layer (PML) absorbing the boundary conditions is used in the program and the corresponding results are discussed. Various types of sources are realized in the program, the results of the waves propagation are given. The test calculations are made for harmonic sources (5 meters - the width of the 2d grid; 100x100 cells; 5GHz frequency of each source), for a horn antenna and for a plane wave propagating in conductors (grid 1 meter wide; 200x200 cells; with conductivity 0.05 siemens/meter, size of the conductor equivalent to 30 grid cells and with the conductor of one cell size and different conductivities: 0.05, 2, 10, 100 siemens/meter) and dielectric (with relative permittivity equal to 4). The program is reliable and has wide possibilities for further development.
Key words: Maxwell equations, Iee algorithm, electromagnetic waves, TM-wave, PML, TFSF.
1. Введение
Уравнения Максвелла являются одними из основопологающих в физике, и поэтому применяются в огромном спектре прикладных задач, таких как сейсморазведка, ускорение
плазмы и фотонные кристаллы. Поэтому, разумеется, разработано множество численных методов, с помощью которых можно решать систему уравнений Максвелла и которые постоянно совершенствуются. Пример одного из таких методов можно найти в [4] (рассматривается система уравнений МГД). В данной работе рассматривается алгоритм Йи (также называемый FDTD method) как метод решения системы уравнений Максвелла, который обладает рядом преимуществ, таких как бездивергентность (в большинстве других методов возникает численная ошибка: уравнение div5 = 0 не выполняется точно, в связи с чем приходится вводить калибровку и решать дополнительно эллиптическую систему, что весьма затратно), простая декомпозиция по данным (сетка очень простая и данные легко хранятся) и, в частности, простота в построении радиолокационных портретов объектов. Большой интерес также представляет задача создания солвера с помощью метода FDTD на произвольных сетках (не обязательно прямоугольных), которая на данный момент еще не решена. Алгоритм в настоящей работе реализован в виде программного модуля для моделирования распространения электромагнитных волн с различными граничными условиями в 2D- и 3D-случаях, в том числе для неоднородной среды. В качестве примера работы алгоритма приводится расчет ряда модельных задач.
Этой тематикой занимаются в ИПМ РАН им. Келдыша (ускорение плазмы), вместе с тем решая систему уравнений МГД [5]. Интерес представляет интеграция метода FDTD в существующий солвер для оптимизации решения некоторых задач. В моделировании фотонных кристаллов также были достигнуты некоторые результаты: в основном программы для моделирования фотонных кристаллов были написаны на языке MatLab, но являются несовершенными и содержат ошибки. Кроме того, сложности для решения представляет обратная задача: какую структуру необходимо создать, чтобы свет проходил в нужных направлениях и в нужных пропорциях. Часто размеры моделируемой области оказываются достаточно большими, и необходимы огромные объемы вычислений, требующие оптимизации вычислений и параллельности программы [4].
Также существует родственный метод FDFD [6], основанный на переходе в частотную область. В дальнейшем имеет смысл совмещение этих двух алгоритмов для оптимизации программы и расширения круга решаемых задач.
2. Численный метод
Система уравнений Максвелла:
дВ -
— = -V х Е - М, ot
dD
= v х Н - J,
at
VD = 0,
У Б = 0,
где Е - напряженность электрического поля, И - электрическая индукция, Н - напряженность магнитного поля, В - магнитная индукция, .] - плотность электрического тока, М - аналогичная плотность магнитного тока. Про нее можно почитать в [3]. Для линейных, изотропных, недиспергирующих материалов:
3 = еео Е,
В =
где е - диэлектрическая проницаемость, £о - диэлектрическая проницаемость вакуума, ц - магнитная проницаемость, цо - магнитная проницаемость вакуума.
Н = аЕ,
М = а*Н,
где а - электрическая проводимость, а* - аналогичная магнитная проводимость. Таким образом,
™ = -1V >Е - 1а-Я,
от ц ц
™ = 1V хЯ - 1аЕ.
от е £
Рис. 1. Сетка метода FDTD
Рис. 2. Поля в одномерном случае
Суть метода Finite-Difference Time-Domain заключается в том, что электрическое и магнитное поля смещены на сетке друг относительно друга так, что каждая компонента
одного поля окружена компонентами другого поля (см. рис. 1). К примеру, в одномерном случае на первом шаге по времени задано поле Е во всех точках. С помощью численных уравнений нашего метода на следующем шаге по времени рассчитывается поле Н во всех ячейках сетки, причем каждая компонента поля Н окружена двумя компонентами поля Е (в двумерном случае, к примеру, таких компонент будет уже 4). На третьем шаге по времени мы уже рассчитываем поле Е, зная поле Н на предыдущем шаге, и т.д. (см. рис. 2).
Далее будет показан метод получения разностных уравнений. Подробный вывод можно найти в [2].
Пространственные и временные производные дискретизируются их конечно-разностными аналогами:
^(гбх, ]5у, кбг, п51) = К+1/2'3'к х ^1/2'3'к + 0[(АЖ)2], ох дх
п+1/2 п-1/2
^ (г5х^5у,к5г,п51) = ^ ~и^к + 0[(А1)2].
Но при такой дискретизации возникает проблема: поле в правой части уравнений для полей будет считаться в момент времени п (производные по координате), а в левой - в моменты п-1/2 и п+1/2. Решим данную проблему, представив поле в точке п в следующем виде:
Ех\п+1/2 + Ех\п-1к/2
тр \п _ х \г,з,к х \ 1,з,к
= 2 . В конечном счете имеем для компонент полей:
р \п+1/2 = г \ г \п-1/2 +
^хЧ,]+1/2,к+1/2 = Га\г,]+1/2,к+1/2^х\1^+1/2,к+1/2 +
+СЬ \ 1,]+1/2,к+1/2(Нг \ ^+1,к+1/2 - Нг\^,к+1/2 - Ну\Ь+1/2,к+1 + Ну\Ь+1/2,к),
р \п+1/2 = Г \ Р \П-1/2 +
^УН-1/2,з+1,к+1/2 = Га\г-1/2,]+1,к+1/2^у\1-1/2^+1,к+1/2 +
+СЬ \ 1-1/2,]+1,к+1/2(Нх \ ™-1/2,]+1,к+1 - Нх\™-1/2,]+1,к - Н*\Хз+1,к+1/2 + ^^\^-1,з+1,к+1/2),
Г \п+1/2 = Г \ Г \п-1/2 +
г-1/2,з+1/2,к+1 = Га\г-1/2,]+1/2,к+1^г\.1-1/2^+1/2,к+1 +
+СЬ \ 1-1/2,]+1/2,к+1(Ну \ "1^+1/2^+1 - Ну\Г1-1,]+1/2,к+1 - Нх\Г1-1/2,]+1,к+1 + Нх\'1-1/2^,к+1),
Нх\™-+1/2^+1,к+1 = ^а\г-1/2,]+1,к+1^х\\^-1/2,]+1,к+1 +
+ П| ( Г \п+1/2 \П+1/2 \П+1/2 \П+1/2 )
+^Ь \ г-1/2,з+1,к+1(^у \ 1-1/2,з+\,к+3/2 ^УЧ-1/2,з+1,к+1/2 ^^1-1/2,з+3/2,к+1 + ^^1-1/2,з+1/2,к+1),
Ну\7++1/2,к+1 = Ва\г,3+1/2,к+1Ну\?Л+1/2,к+1 +
+ П I ( Г \п+1/2 _ Г \п+1/2 _ Г \п+1/2 + Г \п+1/2 )
\ г,з+1/2,к+1 \ 1-1/2,]-1/2,к+1 ^гН-1/2,з+1/2,к+1 ^х\1,з+1/2,к+3/2 + ^х\1,з+1/2,к+1/2),
Н2\7++1,к+1/2 = Ва\г,3+1,к+1/2Нг\™2+1>к+1/2 + 1/2
-3/2,к+1
где (размер ячейки сетки А):
+ П I ( Г \п+1/2 _ Г \П+1/2 _ ^ ¡п+1/2 + Г \П+1/2 )
\ г,з+1,к+1/2(^х \ г,з+з/2,к+1/2 ^х1,з+1/2,к+1/2 ^УН+1/2,з+1,к+1/2 + ^У\1-1/2,з+1,к+1/2),
&1,3,кАА , ( &г,3,кА£
- -зйт; " 1+
Л + А \
V 2£1,3,к ) '
Ъ|г,,?,к (е^л) ^ ( 2£г,.л ) '
а( 2/ЧЛ ) ^ 2)
= I 1 - , ,
2^г/ \
Важными для рассмотрения являются случаи ТМ- и ТЕ-волны. Для них численные уравнения будут иметь следующий вид. ТМ-волна:
р |™+1/2 _ п I р |га-1/2
Е|г-1/2,^+1/2 - °«|г-1/2,^+1/2^г|г-1/2,^+1/2 +
+Съ|г-1/2,,?+1/2(#у ^.7+1/2 - НУ к-^^//™ + - Нх |^_1/2,^+1)'
\п+1
1/2,^+1 — 1 1/2,^+1 ^ж1™-1/2,^+1 +
ТЕ-волна:
+А|г-1/2,.?+1(Е2: |Г-11/22^+1/2 - Е |г-1/2,^+3/2« + 1/2), +А|г,^'+1/2(Ег |Г-+11/22^+1/2 - Е |г-1/2,^+1/2« + 1/2).
ЕУ |Г-1/2,^+1/2 — 1/2,^+1/2 Еу |Г-1/2,^+1/2 + ^Ъ|г-1/2,^+1/2(^г |?-1,^'+1 - Дг |г,^+1П),
^++1/2 — ^¿¿+1/2 Дг ^+1/2 +
+^Ъ|г,^+1/2(Ег |1^++//2 - Е к.7+1/2п + 1/2 + Ег ^Л//^ - Е |г+1/2,^+1« + 1/2).
Как известно, численные методы могут расходиться при неправильном выборе параметров. Кроме того, имеет место численная дисперсия: численное решение может себя вести в принципе не совсем физично: например, сферическая волна может распространяться не одинаково во всех направлениях, а зависеть от направления, в котором распространяется, и тогда получится уже как будто у данной волны форма, похожая чем-то на куб. Подобных неприятностей можно избежать, правильно выбирая шаги по времени и пространству, исходя из следующих соотношений.
Условие сходимости схемы (условие Куранта):
ст I + ^ + ^ < 1 V Ь-х Щ
где с — - скорость света в среде.
Закон численной дисперсии:
(¿4 )2 — (= 2 + (¿4 )2 + (¿4^) )2
Граничные условия
В программе реализовано 3 типа граничных условий: отражающие, периодические и неотражающие.
Самыми простыми в реализации являются отражающие граничные условия. Самый их естественный вид - «идеальный проводник» (perfect electric conductor). Для их реализации достаточно задать на границе значение поля, равное нулю. Другой вариант отражающих граничных условий - имитация столкновения с такой же волной, но идущей в противоположном направлении: задание в граничных ячейках значений поля, равных значению в ячейках у границы. Также можно просто умножить значение поля в необходимых ячейках на очень малую величину (например, 10-10), но это хоть и просто, но не физично, и может быть использовано только в тестовой программе. Последний способ моделирования отражения - задание высокой проводимости (см. а в системе уравнений Максвелла) в необходимых ячейках.
Периодические граничные условия используются в основном для моделирования повторяющихся структур в направлении одной из осей координат (что помогает сделать сетку гораздо меньше и тем самым значительно ускорить вычисления). Суть этих условий - в задании на одной границе значений поля, равных таковым на противоположной границе. Одномерный случай проблем не вызывает. В случаях же большей размерности пространства периодические граничные условия бывает реализовать не всегда просто и даже не всегда возможно: если нам нужно реализовать периодические граничные условия, при этом волна падает под углом на границу области и источник волны пульсирующий, то с реализацией возникают серьезные проблемы. Если же из перечисленных 3 ситуаций выполняются одновременно любые 2 (не все 3), то метод прекрасно работает.
Рис. 3. PML граничные условия
Последний тип граничных условий - неотражающие. В работе используются неотражающие граничные условия, построенные на алгоритме Perfectly Matched Layer [8]. Основной идеей данных граничных условий является определенное задание проводимости на краях сетки в слое толщины более одной ячейки (см. рис. 3). Преимуществом таких поглощающих граничных условий является то, что они верно работают при волнах любой поляризации, частоты и излучаемых любым источником (до 1994 года были также Absorbing Boundary conditions, но они работали только для плоских волн и потому имели крайне ограниченное применение). Данный алгоритм основан на разделении каждой компоненты поля на две ортогональных друг другу (то есть вместо системы из 6 уравнений, записанных ранее для 3-мерного случая, имеем систему из 12 уравнений).
При реализации указанного типа граничных условий есть нюансы, которые необходимо учесть, иначе они будут неверно работать (появится небольшое отражение): проводимость на краях области нужно сделать не постоянной во всех ячейках слоя РМЬ, а возрастающей от центра к границе сетки (например, как ж3). Кроме того, минимизация величины отраженной волны (она может составлять как 10-2, так и 10-5 от величины падающей волны) требует отдельного исследования и тщательного подбора параметров.
В данном случае (как уже было сказано, каждая компонента поля разделяется на две) имеем систему уравнений в 3Б:
д \ д
+ ау\ E%y — Qy (Hzx + Hzy),
д д Sdjj. + az J Exz = — (Hyx + Hyz),
f д \ д
[£dt +a4 Eyz — ~Q~Z(Hxy + Hxz),
д д
+ j Ey% — ~3x(Hzx + Hzy),
д д
+ Ezx — — (Hyx + HyZ),
д д
+°y) EzV — — дд^у (НхУ + Hxz),
( д Л д „ ,
y^Ot +a4 Hxy — — ду( + y),
д Л д
+a*zj Hxz — (Eyx + Eyz),
д д
^ +V*z ) Hyz — — (Exy + Exz),
д д уШ + ay Hyx — дх (Ezx + Ezy),
д д ^Dt + a4 Hzx — — дх(Eyx + Eyz),
д д
+az) Hzy — ду( y + )
Для случая TM- и TE-волн можно довольно просто получить системы уравнений как частные случаи (см. [1]).
Источники излучения
В качестве источников в программе могут быть заданы гармонический источник, гаус-сиан и плоская волна. Задание первых двух не представляет никакой сложности (просто задание поля в определенных ячейках сетки). При моделировании плоской волны все обстоит немного сложнее. Ее можно промоделировать с помощью техники Total Field-Scattered Field [7] (TFSF; см. рис. 4). Основная сложность - необходимость задания источника, излучающего в одном направлении, а не симметрично в разных (в отличие от любого простого источника, заданного значением поля в определенных ячейках). Это делается следующим образом. Расчетная область делится на области Total field (в которой есть как отраженное,
так и поле от источника) и Scattered field (поле только от источника). Граница же между этими двумя областями и является «источником» плоской волны. Для реализации данной идеи задается поле во всех точках границы (например, как гауссиан), а в граничных точках области отраженного поля просто напросто это поле вычитается из имеющихся значений. Поле Н, вычитаемое в области отраженного поля, если поле Е задано как гауссиан g(t), с учетом сдвига по времени и координате между полями:
н = -Л«(
пАу At
t +--- +--
2с 2
) •
Рис. 4. TFSF
3. Результаты
В данном разделе будут показаны результаты некоторых расчетов, проведенных с использованием написанной в процессе исследований программы.
В 2-мерном случае параметры расчетов использовались следующие: частота 5 ГГц, размер расчетной области 5 х 5 м (100 х 100 ячеек сетки; 5 см - размер ячейки). На рис. 5 представлены результаты моделирования работы рупорной антенны. На рис. 6 можно наблюдать распространение излучения от 4 источников, излучающих по гармоническому закону.
Далее на сетке 100 х 500 ячеек с размером ячейки 1 мм представлены примеры моделирования плоской волны (ТЕБЕ) (рис. 7). В области моделирования расположены один за другим диэлектрики с диэлектрической проницаемостью, равной 2 и 4 соответственно, а за ними - проводник с проводимостью 0.1 См/м. Амплитуда плоской волны - 1 В/м; вид импульса показан в одномерном разрезе на рис. 7.
Рис. 5. Рупорная антенна
Рис. 7. Прохождение плоской волны через сложную структуру
4. Заключение
Написана программа, в которой моделируется распространение электромагнитных волн с помощью метода FDTD. В программе можно использовать различные типы источников и граничных условий. Промоделирована работа антенн. Также возможно моделирование неоднородной среды. Данные исследования могут быть применены в широком классе задач. На настоящий момент они применяются в сейсморазведке в двух различных задачах и в фотонных кристаллах. Планируется использование программы в других задачах и областях.
Литература
1. Kane Y. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media // IEEE Transactions on Antennas and Propagation. 1966. V. 14(3). P. 302-307.
2. Taflove A, Hagness S. Computational electrodynamics: the finite-difference time-domain method. 3d edition. London: Artech House, 2005. P. 51-80.
3. Жданов М.С. Электроразведка. М.: Недра, 1986. С. 5-12.
4. Rumpf R.C., Garcia C.R., Berry E.A., Barton J.H. Finite-Difference Frequency-Domain Algorithm for Modeling Electromagnetic Scattering from General Anisotropic Objects // PIERS B. 2014. V. 61. P. 55-67.
5. Holland R., Williams J. Total-field versus scattered-field finite-difference Codes: A Comparative Assessment // IEEE Trans. Nuclear Science, 1983. V. 30(6). P. 4583-4588. DOI: 10.1109/TNS.1983.4333175.
6. Berenger J.-P. A perfectly matched layer for the absorption of electromagnetic waves // J. Computational Physics. 1994. V. 114(2). P. 185-200.
References
1. Kane Y. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Transactions on Antennas and Propagation. 1966. V. 14(3). P. 302-307.
2. Taflove A., Hagness S. Computational electrodynamics: the finite-difference time-domain method. 3d edition. London: Artech House, 2005. P. 51-80.
3. Zhdanov M.S. Electrical Exploration. M.: Nedra. 1986. P. 5-12.
4. Rumpf R.C., Garcia C.R., Berry E.A., Barton J.H. Finite-Difference Frequency-Domain Algorithm for Modeling Electromagnetic Scattering from General Anisotropic Objects. PIERS B. 2014. V. 61. P. 55-67.
5. Holland R., Williams J. Total-field versus scattered-field finite-difference Codes: A Comparative Assessment. IEEE Trans. Nuclear Science, 1983. V. 30(6). P. 4583-4588. DOI: 10.1109/TNS.1983.4333175.
6. Berenger J.-P. A perfectly matched layer for the absorption of electromagnetic waves. J. Computational Physics. 1994. V. 114(2). P. 185-200.
Поступила в редакцию 05.09.2016