CC BY
143
УДК 537.876+004.94 Моделирование распространения электромагнитных волн методом конечных разностей с помощью openEMS
А. А. Шарапова, Д. С. Кулябов
Кафедра прикладной информатики и теории вероятностей Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198
Существующие методы численного моделирования электромагнитного поля в среде обладают, к сожалению, каждый своими недостатками. Авторы поставили себе задачу проанализировать наиболее популярные методы. В качестве модельной задачи авторами рассматривается линза Люнеберга.
В данной работе авторы рассматривают метод конечных разностей во временной области, программное средство openEMS и его применимость к задачам численного моделирования распространения электромагнитных волн в среде на примере сферической линзы Люнеберга.
Благодаря своей простоте и широким возможностям метод конечных разностей во временной области (Finite-Difference Time-Domain method, FDTD) применяется для решения широкого спектра задач. Существует достаточно большое количество программных инструментов, как с открытым исходным кодом, так и проприетарных, позволяющих производить расчёт этим методом. Программный комплекс openEMS является набором функций для MATLAB или Octave, с помощью которых можно произвести расчёт характеристик электромагнитного поля методом EC-FDTD в декартовых или цилиндрических координатах. Программный комплекс openEMS является бесплатным и имеет открытый исходный код. Поддерживаются параллельные технологии вычисления (MPI).
В данной работе на примере моделирования прохождения электромагнитных волн сквозь сферическую линзу Люнеберга показан процесс работы с openEMS, его установка и настройка, а также даны общие сведения о работе метода FDTD и алгоритма Йи. Приведён пример работы алгоритма Йи. Показан способ анализа и визуализации результатов моделирования с помощью программы ParaView. Приведён исходный код скрипта для моделирования.
Исследованы возможности openEMS и метода FDTD при моделировании распространения электромагнитных волн в среде.
Ключевые слова: Метод конечных разностей во временной области, openEMS, линза Люнеберга.
1. Введение
Метод конечных разностей во временной области (далее FDTD) является, возможно, самым простым и популярным методом численной электродинамики. Как и всякий численный метод, он имеет свои достоинства и недостатки. К его достоинствам помимо простоты относится возможность получить результат для широкого спектра длин волн за один расчёт, возможность задать свойства материала в любой точке расчётной сетки, что позволяет задавать анизотропные, дисперсные и нелинейные среды, возможность наблюдать реальное поведение полей во времени и высокая параллельная эффективность. В тоже время FDTD может быть очень ресурсозатратным методом, особенно при моделировании длинных предметов, таких как, например, провода.
В узком смысле под FDTD понимается использование базового алгоритма Йи для численного решения уравнений Максвелла. Разумеется, с момента выхода статьи Йи он был значительно расширен. Сейчас FDTD включает в себя множество возможностей: моделирование сред с дисперсными и нелинейными свойствами, применение различных типов сеток, постропроцессорная обработка результатов и так далее [1,2].
Статья поступила в редакцию 1 января 2016 г.
Работа частично поддержана грантами РФФИ № 14-01-00628, 15-07-08795 и 16-07-00556.
Существует достаточно много программных средств, позволяющих производить расчёт методом FDTD. Одним из них является симулятор openEMS [3]. openEMS не является отдельной программой, а представляет набор функций, которые можно использовать в MATLAB или Octave скриптах. В данной работе будет рассмотрен расчёт электромагнитного поля методом EC-FDTD (Equivalent Circuit FDTD) при помощи openEMS.
Работа имеет следующую структуру. В разделе 2 рассматривается классический алгоритм Йи и его применение для расчёта электромагнитного поля в трёхмерном пространстве, а также приведён стандартный алгоритм расчёта методом FDTD. В разделе 3 описываются некоторые особенности работы с пакетом openEMS. В разделе 4 описывается собственно процесс моделирования сферической линзы Люнеберга с помощью электромагнитного симулятора openEMS.
о
2. Алгоритм Ии
Алгоритм Йи представляет собой схему дискретизации уравнений Максвелла, записанных в дифференциальной форме. Сетки для электрического и магнитного полей смещены по отношению друг к другу на половину шага дискретизации по каждой из пространственных переменных и по времени. В результате узлы, соответствующие компонентам Е, расположены таким образом, что каждый из них окружён четырьмя компонентами Н, и наоборот. Для расчёта значений Е на временном шаге п + 1/2 используются значения Н на шаге п. Аналогичным образом значения Н на шаге п + 1 рассчитываются с использованием значений Е на шаге п + 1/2. Так последовательно рассчитываются все значения полей.
Таким образом, конечно-разностные уравнения позволяют определить электрические и магнитные поля на данном временном шаге на основании известных значений полей на предыдущем, и при заданных начальных условиях вычислительная процедура даёт решение во времени от начала отсчёта с заданным временным шагом.
Ниже описан алгоритм Йи по шагам:
1. Заменить производные в законе Ампера и законе Фарадея конечными разностями.
2. Выразить из полученных дифференциальных уравнений неизвестные значения полей текущего временного шага через уже известные значения полей предыдущего временного шага.
3. Найти значения напряжённости магнитного поля для текущего временного шага.
4. Найти значения напряжённости электрического поля для текущего временного шага.
5. Повторять последние два пункта до тех пор, пока не будут получены значения полей на всём искомом временном промежутке.
2.1. Пример работы алгоритма Йи
Рассмотрим работу метода ЕБТБ поэтапно.
Как уже говорилось, на первом шаге необходимо заменить производные в законах Ампера и Фарадея конечными разностями. Пусть эти уравнения записаны в следующем виде:
^Н
-ат Н - = Ух Е, (1)
В Е
аЕ + е— = Ух Н. (2)
Здесь Н и Е - напряжённости электрического и магнитного полей, ^ и е — диэлектрическая и магнитная проницаемости среды, атН и аЕ — плотности электрического тока и её магнитный аналог.
Обозначим компоненты векторов Н и Е следующим образом: Нх(х,у,х,Ь) = Нх(тАх,пАу ,рАг ,дА1) = Щ [т,п,р], Ну (х,у,г,г) = Ну (тАх,пАу ,рАг ,дА1) = [т,п,р], Нг (х,у,х,Ь) = Нг (тАх,пАу ,рАг = Щ [т,п,р], Ех(х,у,г,г) = Нх(тАх,пАу ,рАг ,дА1) = Е% [т,п,р], Еу (х,у,г,г) = Ну (тАх,пАу ,рАг ,дА1) = [т,п,р], Ег (х,у,г,£) = Нг (тАх ,пАу ,рАг ,дА1) = Е\ [т,п,р].
(3)
(4)
(5)
(6)
(7)
(8)
Здесь х,у,х — координаты узлов в пространстве, а £ — координата по времени. Соответственно Аж, Ау, Аг и А1 — сеточные шаги по соответствующим направлениям.
Как уже было сказано выше, существуют различные типы расчётных сеток, но мы будем пользоваться сеткой Йи как самой распространённой. Расположение узлов в сетке Йи показано на рис. 1.
И2 ЕУ
1 А
1 Лу |
1 1Е2 хТнх Нх 1 1
1 Н2 | Еу
Ех^-/ У-
ЕУ
Ех
Ех
->- У
Рис. 1. Расположение узлов в трёхмерной сетке Йи
Распишем теперь уравнения (1) и (2) покомпонентно:
ЭНХ дЕг дЕу
[Л
—атНг — ^
аЕх + е
ду дх
дНу дЕх дЕг
дх дх
днх = дЕу дЕх
дЬ дх ду '
аЕу + £
(гЕх + £
дЕх днх дНу
т ду дх
дЕу дНх дНг
т дх дх
дЕг дНу дНх
т дх ду
(9)
(10) (11) (12)
(13)
(14)
х
Как можно видеть из этих уравнений, производная по времени каждого компонента выражена через производные по направлениям другого поля.
Теперь можно заменить частные производные конечными разностями. Рассмотрим уравнение (9):
НЧ+ 2
нХ
—
1 1
т,п + 2 -Р +2
+ нХ 2
11
т,п +2-Р +2
7-9+2
Х
2
11
т,п +2-Р +2
-нГ2
11
т,п + 2-Р +2
А,
Е 9
т,п +1,р +2
-Е%
т,п,р + 2
А „
Е9
У
т,п + 2 ,Р + 1
- Ея
У
т,п +- ,р
А „
(15)
Отсюда выразим искомое значение нХ
9+ 2
11
т,п + 2-Р +2
н
9+2
+
11
Р
А
т,п + 2-Р +2
1
2 ц
1+
нХ
2ц
т,п + 2-Р +2
+
Е9
^Ау V У 1
т,п + 2 ,Р + 1
- Е9
у
т,п ^2,Р
А
1+
2 ц
Е9
т,п + 1,р +2
-Е9
т,п,р + 2
. (16)
Точно также следует поступить с уравнениями (10) и (11). После того, как новые значения напряжённости магнитного поля будут найдены, их можно будет использовать для нахождения следующих значений напряжённости электрического поля (уравнения (12)-(14)). Процесс продолжается до тех пор, пока не будут найдены значения Н и Е на всем искомом временном промежутке.
о-т Д*
2.2. Порядок расчёта методом ЕВХЮ
Порядок расчёта методом РБТБ следующий:
1. Задание расчётной области, сетки и граничных условий.
2. В расчётной области задаются материальные тела, составляющие интересующую нас структуру и их свойства.
3. Задание источника излучения.
4. Задание «детектора», отслеживающего изменения полей в нужных точках. Детектор — это не обязательно материальное тело, им может служить просто набор точек.
5. Источник генерирует электромагнитную волну нужного диапазона. Волна падает на тела и рассеивается на них. История распространения волны фиксируется детекторами.
3. Программный пакет моделирования распространения электромагнитных волн openEMS
Программный комплекс openEMS — это электромагнитный симулятор с открытым исходным кодом, использующий для расчёта электромагнитного поля метод EC-FDTD. В отличие от многих других средств оптического моделирования, например, рассмотренных нами в [4], openEMS не является отдельной программой, но набором функций для MATLAB или Octave, написанных на CH—+.
Рассмотрим примеры использования openEMS для решения задачи моделирования распространения электромагнитных волн в различных структурах.
Поскольку для вычислений пакет openEMS использует программные комплексы MATLAB или Octave, необходимо сконфигурировать эти комплексы так, чтобы они могли использовать ресурсы openEMS. Для этого надо добавить путь к библиотекам openEMS (см. [5]):
addpath('~/opt/openEMS/share/openEMS/matlab'); addpath('~/opt/openEMS/share/CSXCAD/matlab');
В случае MATLAB следует добавить строчки в файл startup.m, а в случае Octave — в файл .octaverc.
Визуализировать результаты симуляции можно с помощью программы ParaView. Её можно установить из стандартных репозиториев. В ParaView необходимо открыть папку, в которой сохранен результат работы скрипта, и выбрать .vtr файл (на самом деле их несколько, но в обозревателе они представлены как дерево). Затем нажать кнопку Apply под списком Properties. Чтобы визуализировать амплитуду электромагнитной волны (далее для краткости ЭМВ), нужно в выпадающем списке Coloring выбрать E-field (по умолчанию там стоит SolidColor). Чтобы просмотреть анимацию, нужно нажать кнопку Play на панели инструментов. Также ParaView позволяет анализировать результаты при помощи графиков. Например, зависимость амплитуды ЭМВ от времени в определённой точке расчётной сетки или в зависимости от расстояния [6].
4. Моделирование распространения электромагнитных волн через сферическую линзу Люнеберга
Теперь рассмотрим задачу моделирования сферической линзы Люнеберга в орепЕМБ. Ранее мы уже рассматривали некоторые теоретические аспекты проектирования сферической линзы Люнеберга [7], а также моделирование такой линзы из кубиков, но в рамках геометрической оптики [8,9]. Как известно, относительная диэлектрическая проницаемость в линзе Люнеберга зависит от расстояния до центра линзы и имеет вид
Свойство линзы таково, что она преобразовывает сферический фронт волны, расходящийся из точки на её поверхности, в плоский фронт [10-12].
Пусть линза имеет диаметр 10 Л, где Л — длина волны. Таким образом, чтобы рассчитать диаметр, нужно знать длину волны, но в орепЕМБ можно задать
3.1. Конфигурирование openEMS
3.2. Визуализация с помощью ParaView
4.1. Постановка задачи
(17)
только частоту излучения. Для нахождения оставшегося параметра воспользуемся формулой
А = -. (18)
При задании параметров излучения в орепЕМЯ необходимо иметь ввиду следующее ограничение:
тах (Ах,у,г) < . (19)
Зададим точечный источник излучения на поверхности линзы. Пусть излучение имеет частоту 10 МГц и длина волны равна 30 м (радиоволны).
Значения электрического поля зафиксированы в плоскости .
4.2. Анализ скрипта openEMS
Ранее мы уже рассматривали синтаксис openEMS на примере более простой модели, волновода с двумя металлическими стенками [13]. Несмотря на то, что линза Люнеберга является более сложной структурой, схема описания модели остаётся той же. Рассмотрим скрипт построчно. Задание пространства FDTD c 300 временными шагами:
FDTD = InitFDTD('NrTS',300, 'EndCriteria',0,'0verSampling',50); Задание частоты излучения (10 МГц): freq = 10e6;
Задание скорости света (м/с) для расчёта длины волны:
c0 = 299792458;
Расчёт длины волны (м):
lambda = c0/freq;
Задание излучения:
FDTD = SetSinusExcite(FDTD,freq);
Задание граничных условий (идеально согласованные слои — Perfectly Matched Layer, PML):
BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; FDTD = SetBoundaryCond(FDTD,BC); CSX = InitCSX();
Задание размеров расчётной сетки. В данном случае используется прямоугольная система координат: start_mesh=-200; end_mesh=200; unit = 5;
mesh.x = SmoothMeshLines([start_mesh end_mesh] mesh.y = SmoothMeshLines([start_mesh end_mesh] mesh.z = SmoothMeshLines([start_mesh end_mesh] CSX = DefineRectGrid(CSX, 5, mesh);
Задание параметров линзы:
CSX = AddMaterial( CSX, 'sphere_material');
Сначала диэлектрическая проницаемость материала линзы устанавливается равной единице:
CSX = SetMaterialProperty(CSX,'sphere_material', 'Epsilon', 1);
Здесь устанавливается настоящее распределение диэлектрической проницаемости в линзе. К сожалению, радиус линзы нужно задавать в твёрдых числах. Здесь rho — расстояние до оси :
CSX = SetMaterialWeight(CSX, 'sphere_material', 'Epsilon', ['2-(rho*rho)/22470']); Расчёт радиуса сферы: sphere_radius = lambda/2*10;
Добавление примитива сферы с заданным ранее материалом и найденным радиусом:
CSX = AddSphere(CSX,'sphere_material',2,[0 0 0],sphere_radius);
, unit); , unit); , unit);
Добавление «точечного источника», т.е. очень маленькой площадки, с которой распространяется электромагнитная волна:
CSX = AddExcitation(CSX,'excitation',1,[0 20 0]);
CSX = AddBox(CSX,'excitation',0,[-unit*2 -unit*2 -sphere_radius-unit], [unit*2 unit*2 -sphere_radius]);
Задание плоскости, в которой наблюдается распространение электромагнитной волны:
CSX = AddDump(CSX,'Et','DumpType',0,'DumpMode',0);
CSX = AddBox(CSX,'Et',0,[start_mesh 0 start_mesh],[end_mesh 0 end_mesh]); Запись результатов в xml-файл: Sim_Path = 'point_source'; Sim_CSX = 'point_source.xml';
[status, message, messageid] = rmdir(Sim_Path,'s'); [status, message, messageid] = mkdir(Sim_Path); WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); CSXGeomPlot([Sim_Path '/' Sim_CSX],'--RenderDiscMaterial' );
RunOpenEMS(Sim_Path, Sim_CSX);
Результатом расчётов, кроме численных значений, является также визуализация исследуемого процесса, фрагмент которой можно видеть на рис. 2.
Рис. 2. Прохождение электромагнитной волны через сферическую линзу
Люнеберга
5. Заключение
В данной работе даны сведения о работе метода FDTD: описание классического алгоритма Йи, а также общий порядок расчёта методом FDTD. Осуществлено моделирование прохождения электромагнитной волны через сферическую линзу Люнеберга с помощью программы openEMS. Был написан скрипт, решающий данную задачу. Показана применимость openEMS в моделировании градиентной оптики.
Литература
1. Taflove A., Hagness S. C. Computational Electrodynamics: The Finite-difference Time-domain Method. — Artech House, 2005. — ISBN 978-1-58053-832-9.
2. Schneider J. B. Understanding the Finite-Difference Time-Domain Method. — School of Electrical Engineering and Computer Science Washington State University, 2010. — www.eecs.wsu.edu/schneidj/ufdtd.
3. Liebig T. openEMS - Open Electromagnetic Field Solver. — http://openEMS. de.
4. Шарапова А. А. Обзор бесплатного программного обеспечения для моделирования линз и оптических систем // Программа 57-й научной конференции МФТИ. — МФТИ, 2015.
5. openEMS Online Manual. — http://openems.de/index.php/Introduction.
6. Moreland K. The ParaView Tutorial. — 2015.
7. Шарапова А. А. Нахождение распределения коэффициента преломления для линзы Люнеберга произвольной конфигурации // Информационно-телекоммуникационные технологии и математическое моделирование высокотехнологичных систем: материалы Всероссийской конференции с международным участием. Москва, РУДН, 20-24 апреля 2015 г. — М.: РУДН, 2015. — С. 323-325.
8. Шарапова А. А. Математическое моделирование линзы Люнеберга из кубиков // Труды 55-й научной конференции МФТИ. — МФТИ, 2012. — С. 164.
9. Шарапова А. А. Математическое моделирование сферической линзы Люнеберга из кубиков // Информационно-телекоммуникационные технологии и математическое моде- лирование высокотехнологичных систем: материалы Всероссийской конференции с международным участием. Москва, РУДН, 22-26 апреля 2013 г. — РУДН, 2013. — С. 221-222.
10. Устройства СВЧ и антенны / Д. И. Воскресенский, В. Л. Гостюхин, В. М. Максимов, Л. И. Пономарев; под ред. Д. И. Воскресенский. — Радиотехника, 2006.
11. Greenwood A. D., Jin J.-M. A Field Picture of Wave Propagation in Inhomoge-neous Dielectric Lenses // IEEE Antennas and Propagation Magazine. — 1999.
12. Modeling and simulations of Luneburg lens antennas for communication purposes / S. R. Baev, S. M. Gechev, B. N. Hadjistamov, P. I. Dankov // 16th Telecommunications forum FOR, Serbia, Belgrad. — 2008. — Pp. 488-491. — http://2008.telfor.rs/files/radovi/07_04.pdf.
13. Шарапова А. А. Моделирование распространения электромагнитных волн с помощью openEMS // 58-я научная конференция МФТИ с международным участием. Москва, МФТИ, 23-28 ноября 2015 г. — МФТИ, 2015. — С. 68. — http://conf58.mipt.ru/static/reports_pdf/274.pdf.
UDC 537.876+004.94
Simulation of Wave Propagation with openEMS
A. A. Sharapova, D. S. Kulyabov
Department of Applied Probability and Informatics Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., Moscow, Russian Federation, 117198
Unfortunately, all the existing methods used for modeling computational electrodynamics have their weaknesses. The authors' goal is to analyze the most popular methods. We use spherical Luneburg lens as an illustration.
In this paper authors review the Finite-difference time-domain (FDTD for short), electromagnetic field solver openEMS and its applicability for simulation of wave propagation through medium with the spherical Luneburg lens as an example.
Thanks to it's simplicity and broad capabilities, the FDTD method is widely used in various fields. There are quite a lot of simulation tools that implement FDTD, both open-source and proprietary. OpenEMS is an extension for MATLAB and Octave for solving electromagnetic field using the EC-FDTD method. It supports Cartesian and cylindrical coordinates. OpenEMS is free and open-source. It also supports multi-threading, SIMD (SSE) and MPI.
In this paper we simulate wave propagation through spherical Luneburg lens via openEMS thus showing the capabilities of the tool, such as a method to simulate GRIN-optics. We also show how to install and configure it, and how to visualize and analyze the results using ParaView, the application for scientific visualization. Source code of the simulation is presented with appropriate commentaries. In this paper we also briefly review the classic FDTD method and Yee algorithm.
We examined the capabilities of openEMS and FDTD method to simulate wave propagation in a medium.
Key words and phrases: Finite difference time domain method, openEMS, Luneburg lens.
References
1. A. Taflove, S. C. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method, Artech House, 2005.
2. J. B. Schneider, Understanding the Finite-Difference Time-Domain Method, School of Electrical Engineering and Computer Science Washington State University, 2010.
URL www.eecs.wsu.edu/schneidj/ufdtd
3. T. Liebig, openEMS - Open Electromagnetic Field Solver. URL http://openEMS.de
4. A. A. Sharapova, Review of Free Software for Modeling of Lenses and Optical Systems, in: 57th MFTI Scientific Conference, MFTI, 2015, in Russian.
5. openEMS Online Manual.
URL http://openems.de/index.php/Introduction
6. K. Moreland, The ParaView Tutorial (2015).
7. A. A. Sharapova, Calculation of Distribution the Refractive Index in a Luneburg Lens of an Arbitrary Configuration, in: Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems, Moscow, PFUR, 20-24 April 2015, PFUR, 2015, in Russian.
8. A. A. Sharapova, Numerical Modeling of a Luneburg Lens Formed from Cubes, in: 55th MFTI Scientific Conference, MFTI, 2012, p. 164, in Russian.
9. A. A. Sharapova, Numerical modeling of a spherical luneburg lens formed from cubes, in: Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems, PFUR, 2013, in Russian.
10. D. I. Voskresenskij, V. L. Gostyuhin, V. M. Maksimov, L. I. Ponomarev, Microwave Devices and Antennas, Radio Engineering, 2006, in Russian.
11. A. D. Greenwood, J.-M. Jin, A Field Picture of Wave Propagation in Inhomogeneous Dielectric Lenses, IEEE Antennas and Propagation Magazine.
12. S. R. Baev, S. M. Gechev, B. N. Hadjistamov, P. I. Dankov, Modeling and Simulations of Luneburg Lens Antennas for Communication Purposes, in: 16th Telecommunications forum FOR, Serbia, Belgrad, 2008, pp. 488-491.
URL http://2008.telfor.rs/files/radovi/07_04.pdf
13. A. A. Sharapova, Simulation of Wave Propagation with openEMS, in: 58th MFTI Scientific Conference. Moscow, MFTI, 23-28 November 2015, MFTI, 2016, p. 68, in Russian.
