Диссертация (Компьютерная реализация геометрических методов в максвелловской оптике), страница 13

Порт необязательно должен быть привязан к примитиву.Для расчета значений напряжения, силы тока и некоторых другиххарактеристик на порте используется функция calcPort. СинтаксисcalcPort имеет вид :port = calcPort(port, SimDir, f, varargin)154port — имя для возвращаемого объекта, SimDir — директория, в которой лежат файлы симуляции, f — вектор частот для дискретногопреобразования Фурье, varargin — дополнительные параметры, зависящие от вида порта.Нами были рассмотрены все основные функции и концепции, необходимые для моделирования в openEMS.

В следующем разделе намибудет рассмотрена визуализация и анализ результатов.Визуализация с помощью ParaViewВизуализировать результаты симуляции можно с помощью программы ParaView [33]. Ее можно установить из стандартных репозиториев.В ParaView необходимо открыть папку, в которой сохранен результатработы скрипта, и выбрать .vtr файл (на самом деле их несколько, нов обозревателе они представляются как дерево). Затем нажать кнопкуApply под списком Properties.

Чтобы визуализировать амплитуду ЭМВ,нужно в выпадающем списке Coloring выбрать E-field (по умолчаниютам стоит SolidColor). Чтобы просмотреть анимацию, нужно нажатькнопку Play на панели инструментов. Также ParaView позволяет анализировать результаты при помощи графиков. Например, зависимостьамплитуды ЭМВ от времени в определенной точке расчетной сетки илив зависимости от расстояния.155C.2.3.Моделирование распространения электромагнитныхволн в сферической линзе Люнеберга с помощью пакетаopenEMSТеперь рассмотрим задачу моделирования сферической линзы Люнеберга [69] в openEMS [151 ; 185].Как известно, относительная диэлектрическая проницаемость в линзе Люнеберга зависит от расстояния до центра линзы и имеет видεr (r) = 2 − r 2R(0 6 r 6 R).Свойство линзы таково, что она преобразовывает сферическийфронт волны, расходящийся из точки на её поверхности, в плоскийфронт.Значения электрического поля зафиксированы в плоскости XZ.Анализ скрипта моделирования сферической линзы ЛюнебергаРассмотрим упрощённый скрипт для моделирования.Задание пространства FDTD c 300 временными шагами :FDTD = InitFDTD('NrTS',300,,→'EndCriteria',0,'OverSampling',50) ;Задание частоты излучения (10 МГц) :freq = 10e6;Задание скорости света (м/с) для расчёта длины волны :156(a) Сечения фронта волны в линзеЛюнеберга, t=50 с(b) Сечения фронта волны в линзеЛюнеберга, t=100 с(c) Сечения фронта волны в линзеЛюнеберга, t=150 с(d) Сечения фронта волны в линзеЛюнеберга, t=200 сРис.

C.1. Сечения фронта волны в линзе Люнеберга157c0 = 299792458;Расчёт длины волны (м) :lambda = c0/freq;Задание излучения :FDTD = SetSinusExcite(FDTD,freq);Задание граничных условий (PML — идеально согласованные слоиили perfectly matched layer) :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], unit) ;mesh.y= SmoothMeshLines([start_mesh end_mesh], unit) ;mesh.z= SmoothMeshLines([start_mesh end_mesh], unit) ;CSX = DefineRectGrid(CSX, 5, mesh) ;Задание параметров линзы :CSX = AddMaterial( CSX, 'sphere_material') ;Сначала диэлектрическая проницаемость материала линзы устанавливается равной единице :158CSX = SetMaterialProperty(CSX,'sphere_material', 'Epsilon',1) ;Здесь устанавливается настоящее распределение диэлектрической про,→ницаемости в линзе.

К сожалению, радиус линзы нужно задавать втвердых числах. Здесь rho — расстояние до оси z :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) ;Добавление «точечного источника», т.е. очень маленькой площадки,с которой распространяется электромагнитная волна :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-файл :159Sim_Path = 'tmp' ;Sim_CSX = 'tmp.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) ;Результат расчёта линзы Люнеберга показан на рис.

C.1.C.3.Численное решение уравнения эйконала мето-дом характеристикУравнение эйконала [10 ; 116] описывает распространение волн вслучае аппроксимации волнового уравнения с помощью квазиклассического приближения. Уравнение эйконала связывает волновую оптику сгеометрической оптикой.C.3.1.Характеристики уравнений в частных производныхРассмотрим уравнениеa1 (x, y)∂u(x, y)∂u(x, y)+ a2 (x, y)= f (x, y).∂x∂y160Данное уравнение эквивалентно утверждению о том, что векторное поле с компонентами (a1 (x , y) , a2 (x , y) , f (x , y)) является касательным кповерхности z = u(x , y), которая имеет нормальный вектор с компонентами (ux , uy , −1).

Уравнение характеристик записывается в следующем виде :dx= a1 (x, y),dtC.3.2.dy= a2 (x, y),dtdu(x, y)= f (x, y).dtУравнение эйконалаРассмотрим уравнение эйконала с граничным условием|∇u(x)|2 = b2 (x), x ∈ Rn ,u(x) = ϕ(x),x ∈ Γ ⊂ Rn .В двумерном случае2 2∂u(x, y)∂u(x, y)+= b2 (x, y), (x, y) ∈ R2 ,∂x∂yu(x, y) = ϕ(x, y),(x, y) ∈ Γ ⊂ R2 .Для сведения данного уравнения к системе ОДУ, воспользуемся методом характеристик. Для этого следует преобразовать уравнение эйконала к виду (C.3.1) Для этого сделаем следующую замену переменных :p1 =∂u,∂xp2 =∂u.∂y161В этом случае само уравнение преобразуется к виду :|p|2 = p21 + p22 = b2 (x, y).Проделаем далее ряд преобразований.∂ 2∂p1∂p2∂b(p1 + p22 ) = 2p1+ 2p2= 2b ,∂x∂x∂x∂x∂ 2∂p1∂p2∂b(p1 + p22 ) = 2p1+ 2p2= 2b .∂y∂y∂y∂yПолучим систему уравнений∂p∂b∂p1∂p2∂bp,=b ,p1+ p2=b ,∂x∂x∂x∂x ⇒ ∂x ∂p1∂p2∂b∂p∂bp1+ p2=b ,p,=b ,∂y∂y∂y∂y∂yТак как∂p1∂ 2 u(x, y) ∂ 2 u(x, y) ∂p2===,∂y∂y∂x∂x∂y∂xто∂p1∂p2=.∂y∂xИспользуя это равенство преобразуем наши выражения ∂p∂p1 ∂p2∂p1 ∂p1=,=,=∂x∂x ∂x∂x ∂y ∂p∂p1 ∂p2∂p2 ∂p2=,=,=∂y∂y ∂y∂x ∂y∂p1= ∇p1 ,∂x∂p2= ∇p2 .∂x162В результате :∂p∂bp,=b ,∂x∂x∂p∂bp,=b ,∂y∂y⇒∂b,∂x∂b(p, ∇p2 ) = b .∂y(p, ∇p1 ) = bТаким образом цель достигнута — уравнение (C.3.2) преобразованоу двум уравнениям вида (C.3.1).∂p1∂p1∂bp1 ∂p1 p2 ∂p1 1 ∂bp1+ p2=b ,+ 2=, 2∂x∂y∂xb ∂xb ∂y b ∂x⇒∂p2∂p2∂bp ∂pp ∂p1 ∂b 21 2 + 22 2 =p1+ p2=b ,.∂x∂y∂yb ∂xb ∂y b ∂yДля каждого уравнения выпишем характеристикиp1 ∂p1 p2 ∂p11 ∂b+ 2=,2b ∂xb ∂yb ∂xdx p1= 2,dtbdyp2= 2dtbdp11 ∂b=.dtb ∂x163p1 ∂p2 p2 ∂p21 ∂b+=,b2 ∂xb2 ∂yb ∂ydx p1= 2,dtbdyp2= 2,dtbdp21 ∂b=.dtb ∂yИтак мы получили систему ОДУ из четырех уравнений для четырехфункций : x(t), y(t), p1 (t), p2 (t) :dx p1= 2,dtbdyp2= 2,dtbdp11 ∂b=,dtb∂xdp21 ∂b=.dtb ∂yНачальные условияx(t)|t=0 = x0 ,y(t)|t=0 = y0p1 (t)|t=0 = c1 b(x0 , y0 )p2 (t)|t=0 = c2 b(x0 , y0 )На выбор констант c1 и c2 накладывается условие c21 + c22 = 1.Установим связь между параметром t и функцией u(x , y).

Так как :du ∂u dx ∂u dy=+,dt∂x dt ∂y dtи∂u ∂u,∂x ∂y= ∇u = p,164тоdudx= ∇u=dtdtdxdxdyp1 p1 p2 p 2|p|2p,= p1+ p2= 2 + 2 = 2 .dtdtdtbbbВ силу того, что |p|2 = b2 (x , y), получаем :du |p|2b2du= 2 = 2 =1 ⇒= 1.dtbbdtРешением уравнения ut = 1 является функция u(x , y) = t+const из чегоследует, что параметр t имеет физический смысл времени прохождениясигнала от точки (x0 , y0 ) до точки (x , y)C.3.3.Метод FSMСледует учитывать два ключевых момента при численном решенииуравнения эйконала :– Введение устойчивой и точной дискретизации уравнения в частныхпроизводных.– В связи с нелинейностью задачи приходится решать большую систему нелинейных уравнений, поэтому необходим эффективный методее решения.Выделяют два подхода к численному решению уравнения эйконала :– Преобразование к системе обыкновенных дифференциальных уравнений (системе уравнений Гамильтона), а затем применение одного165из многочисленных методов численного решения таких уравнений.– Подход к задаче как с стационарной краевой задаче.

Разработкаэффективного численного алгоритма решения системы нелинейныхуравнений, получившихся при дискретизации. К этому типу методов относится, например, метод быстрой прогонки (fast marchingmethod).Метод fast sweeping method, FSM (метод быстрого подметания/уборки) был предложен достаточно недавно [26 ; 112]. Основная идея методазаключается в использовании противопотоковой разностной схемы Годунова и итерационной схемы Гаусса-Зейделя с переменным порядкомпрохода узлов сетки.FSM прост в реализации и требует конечного числа итераций. Сложность алгоритма O(N ) для N точек сетки. Число итераций не зависитот числа узлов сетки (от размера сетки).

FSM метод можно распространить на общий случай уравнения Гамильтона-Якоби [133].Рассмотрим решение двумерного уравнений эйконала [117 ; 135 ;156].FSM для двумерного уравнения эйконалаРассмотрим уравнение∂u∂x2+∂u∂y2= b2 (x, y) > 0166yJyjy1Ox1xixIРис. C.2. Различные точки сетки разбиения области интегрированияс граничным условием u(x , y) = 0, (x , y) ∈ Γ ⊂ R2 . Опишем для этогоуравнения алгоритм F SM .Зададим на области интегрирования прямоугольную сетку (см.рис. C.2). Введем I точек разбиения x1 < x2 < x3 < . .