Различия

Здесь показаны различия между выбранной ревизией и текущей версией данной страницы.

ru:tutorial [2012/03/01 12:22]
deinega [Прохождение через плоскопараллельную пластинку]
ru:tutorial [2014/01/11 06:04] (текущий)
deinega [Наклонное падение]
Строка 33: Строка 33:
 ======Геометрические тела====== ======Геометрические тела======
  
-Ниже мы вкратце перечислим классы, отвечающие некоторым используемым геометрическим телам. Со всем списком можно ознакомиться в справке для пользователя. +Как создавать геометрические тела описано [[/ru/ivutils#geometricheskie_tela|здесь]].
- +
-Класс вектор Vector_3. В его конструкторе задаются три координаты x, y и z. Вектора можно складывать, вычитать, умножать на число, умножать скалярно (с помощью оператора *) и векторно (с помощью оператора %). Например, +
- +
-<code cpp> +
-  Vector_3 v1(1,0,0),v2(0,1,0); +
-  Vector_3 v0=a%b-Vector_3(0,0,1); // v0 равен Vector_3(0,0,0) +
-</code> +
- +
-Если в конструкторе вектора задается меньше трех координат, то незаданные координаты полагаются равными последней заданной. Например, Vector_3(1)=Vector_3(1,1,1). Также Vector_3()=Vector_3(0,0,0), так как первая координата по умолчанию равна нулю. +
- +
-Векторам с целочисленными координатами отвечает класс iVector_3. +
-   +
-Прямоугольный параллелепипед, ребра которого параллельны координатным осям, задается путем указания левого нижнего и правого верхнего векторов. +
- +
-<code cpp> +
-  Box b(Vector_3(),Vector_3(1)); // единичный куб b +
-</code> +
- +
-Сфера задается путем указания ее радиуса и центра. +
- +
-<code cpp> +
-  Sphere s(1,Vector_3()); // сфера s радиуса 1, центр которой находится в начале координат +
-</code> +
- +
-Также имеются многогранник, цилиндр и пирамида с основанием в виде произвольного многоугольника. Тела можно растягивать, поворачивать, брать их пересечения и объединения. Гладкие тела, например конус, можно конструировать как многогранники с большим количеством граней.+
  
 Если тела выходят за границы вычислительного объема, то рассматривается только та их часть, которая в него попадает. Если тела выходят за границы вычислительного объема, то рассматривается только та их часть, которая в него попадает.
Строка 185: Строка 160:
 Форма Беренгера выглядит так: Форма Беренгера выглядит так:
  
-<tex>E(t) = E_0 (t-t_0) \exp \left[ \left(\frac{t-t_0}{t_w} \right)^2 \right] </tex>, <tex>\left|t-t_0\right| \le T</tex>+<tex>E(t) = E_0 (t-t_0) \exp \left[ -\left(\frac{t-t_0}{t_w} \right)^2 \right] </tex>, <tex>\left|t-t_0\right| \le T</tex>
  
 <tex>E(t) = 0</tex>, <tex> \left|t-t_0\right| > T </tex> <tex>E(t) = 0</tex>, <tex> \left|t-t_0\right| > T </tex>
  
-В нашем примере <tex>E_0=1</tex>,<tex>t_0=25</tex>,<tex>t_w=5</tex>,<tex>T=25</tex>+В нашем примере <tex>E_0=1</tex>,<tex>t_0=25</tex>,<tex>t_w=5</tex>,<tex>T=125</tex>
  
 =====Детекторы===== =====Детекторы=====
Строка 452: Строка 427:
 Значение "c" и приведут к тому, что будет выполнено построение геометрии и расчет. Это полезно, если вы хотите отложить анализ больших двоичных файлов на потом (например, когда задача будет перезапущена на другом числе процессоров). Значение "c" и приведут к тому, что будет выполнено построение геометрии и расчет. Это полезно, если вы хотите отложить анализ больших двоичных файлов на потом (например, когда задача будет перезапущена на другом числе процессоров).
  
-Значение "a" означает выплонение только фазы анализа.+Значение "a" означает выполнение только фазы анализа.
  
-======Прохождение через плоскопараллельную пластинку======+======Спектры прохождения и отражения для периодических структур======
  
-С помощью FDTD можно рассчитывать спектры прохождения и отражения от различных структур. +С помощью FDTD можно рассчитывать спектры прохождения и отражения для периодических структур (фотонные кристаллы, массивы антенн и т.д.) 
-Здесь мы проиллюстрируем как это делать на примере бесконечной плоскопараллельной пластинки.+Для их расчета рассчитывается распространение ограниченного по времени импульса через структуру. 
 +Периодичность моделируется путем задания периодических граничных условий в нужных направлениях. 
 + 
 +Рассмотрим элементарную ячейку фотонно-кристаллической пластинки, состоящей из сфер упакованных в квадратную решетку. 
 +Такая ячейка ограничена периодическими границами 1 и 2. 
 +Граница Total Field / Scattered Field (TF/SF) 3 генерирует падаюший импульс (область Total Field заштрихована). 
 +Для непериодического направления используются поглощающие граничные уловия PML. 
 +Они поглощают отраженную и прошедшую волны, моделируя их уход на бесконечность. 
 +Численный эксперимент ведется до тех пор пока излучение не выйдет из структуры. 
 +Прошедшая и отраженная волны записываются во время эксперимента, преобразуются в частотное представление и нормализуются на спетр падающей волны, что дает спектры прохождения и отражения. 
 + 
 +{{tutorial_periodic:geometry_norm.png?250}}\\ 
 + 
 +В случае нормального падения периодические граничные условия выглядят так 
 + 
 +<tex> 
 +{\bf F}({\bf x},t)={\bf F}\left({\bf x} \pm {\bf a},t\right), 
 +</tex> 
 + 
 +где <tex>\bf F</tex> это электрическое или магнитное поле (<tex>\bf E</tex> или <tex>\bf H</tex>), <tex>\bf a</tex> это вектор трансляции для периодической структуры, <tex>\bf x</tex> и <tex>t</tex> пространственная и временная координаты. 
 + 
 +В случае наклонного падения периодические граничные условия содержат временной сдвиг.  
 +В 2D они принимают форму 
 + 
 +<tex> 
 +{\bf F}({\bf x},t)={\bf F}\left({\bf x}\pm{\bf a},t\pm a\sin{\theta}/c\right), 
 +</tex> 
 + 
 +где <tex>\theta</tex> это угол падения, <tex>c</tex> - скорость света в среде, откуда падает свет. 
 +Смысл этого выражения легко понять из следующих соображений. 
 +Падающий импульс приходит на границу 1 раньше чем на границу 2, поэтому значения полей на этих границах "сдвинуты во времени" относительно друг друга. 
 +Использование периодических граничных условий требует знания поля в прошлом на границе 1 (для применения на границе 2) и поля в будущем на границе 2 (для применения на границе 1). 
 +Поле в прошлом могут быть получено, если ведется запись распространения волны. 
 +Получение поля в будущем проблематично, поскольку оно неизвестно заранее. 
 +Для решения этой проблемы можно использовать итерационный метод для моделирования наклонного падения на периодическую структуру. 
 + 
 +Ниже мы рассматрим моделирование нормального падения на одномерном примере, далее мы объясним принцир итерационного метода и покажем как он наботает на примере фотонно-кристаллической пластинки из сфер и антиотражающего текстурированного покрытия. 
 + 
 + 
 +=====Бесконечная пластинка===== 
 + 
 +Ниже мы рассмотрим как рассчитать спектры прохождения и отражения для бесконечной пластинки (которая может рассматриваться как периодическая с периодом нуль) в случае нормального падения.
  
 Создадим объект этого класса uiExperiment:  Создадим объект этого класса uiExperiment: 
Строка 535: Строка 551:
   M.AddPolarizability(emLorentz(.5*(2*M_PI), .1*(2*M_PI), 2*1e-5));   M.AddPolarizability(emLorentz(.5*(2*M_PI), .1*(2*M_PI), 2*1e-5));
  
-  task.AddObject(M,GetPolyhedronPlate(Vector_3(0,0,1),Vector_3(0,0,1),1));+  task.AddObject(M,GetPlate(Vector_3(0,0,1),Vector_3(0,0,1),1));
  
-Пластинка задается функцией GetPolyhedronPlate, первый аргумент которой задает направление вдоль пластинки (направление z), второй точку, через которую проходит одна из сторон пластинки, а третий - толщину пластинки (она равна 1).+Пластинка задается функцией GetPlate, первый аргумент которой задает направление вдоль пластинки (направление z), второй точку, через которую проходит одна из сторон пластинки, а третий - толщину пластинки (она равна 1).
 Пластинка сделана из дисперсной среды, описываемой двумя членами Лоренца. Пластинка сделана из дисперсной среды, описываемой двумя членами Лоренца.
  
Строка 554: Строка 570:
 Поглощение может быть получено как 1 - R - T. Поглощение может быть получено как 1 - R - T.
  
-Для задания полубесконечной подложки мы можем воспользоваться функцией GetPolyhedronPlane. В отличие от пластинкиу подложки нет второй замыкающей поверхности, и она бесконечно продолжается вдоль оси z. В модели FDTD она погружается в PMLчто моделирует ее бесконечность.+Расстояние между теламиграницей TF/SFPML и детекторами должно быть не менее 2 сеточных шагов. В нашем примере мы придерживались большего расстояния (5-20 сеточных шагов)поскольку иногда это может улучшить результат.
  
-Материалом подложки будет кремний. Для получения соответствующего объекта emMedium можно воспользоваться функцией getSi, у которой имеется аргумент, по умолчанию равный единице. Этот аргумент задает, чему равна единица длины FDTD в микронах. Это необходимо для правильной калибровки коэффициентов, участвующих в расчете диэлектрической проницаемости. В нашем случае единица FDTD будет равна 0.1 микрона.+=====Кремниевая подложка=====
  
-  task.AddObject(getSi(.1),GetPolyhedronPlane(Vector_3(0,0,1), Vector_3(0,0,1)));+Ниже мы рассчитываем отражение от полубесконечной кремниевой подложки в случае нормального падения. Мы используем уже имеющийся код из предыдущей подсекции (задание граничных условий, области Total Field и т. д.) 
 + 
 +Для задания полубесконечной подложки мы можем воспользоваться функцией GetHalfSpace. В отличие от пластинки, у подложки нет второй замыкающей поверхности, и она бесконечно продолжается вдоль оси z. В модели FDTD она погружается в PML, что моделирует ее бесконечность. 
 + 
 +Материалом подложки будет кремний. Для получения соответствующего объекта emMedium можно воспользоваться функцией getSi, у которой имеется аргумент, по умолчанию равный единице. Этот аргумент задает, чему равна единица длины FDTD в микронах. Это необходимо для правильной калибровки коэффициентов, участвующих в расчете диэлектрической проницаемости (см. [[tutorial?&#sredy|как задавать среды в EMTL]]). В нашем случае единица FDTD будет равна 0.1 микрона. 
 + 
 +  task.AddObject(getSi(.1),GetHalfSpace(Vector_3(0,0,1), Vector_3(0,0,1)));
      
 Параметры, используемые для апроксимации диэлектрической проницаемости кремния функцией getSi, были получены с помощью программы на MatLab, которую вы можете найти в разделе [[fitting]]. Параметры, используемые для апроксимации диэлектрической проницаемости кремния функцией getSi, были получены с помощью программы на MatLab, которую вы можете найти в разделе [[fitting]].
Строка 573: Строка 595:
 {{tutorial_plate:r_si.png}} {{tutorial_plate:r_si.png}}
  
-Как видно, кремний достаточно сильно отражает в оптическом диапазоне. Для уменьшения отражения могут использоваться текстурированные антиотражающие покрытия.+Как видно, кремний достаточно сильно отражает в оптическом диапазоне. Для уменьшения отражения могут использоваться [[ar|антиотражающие текстурированные покрытия]]. 
 +=====Наклонное падение===== 
 + 
 +Итерационный метод для моделирования наклонного падения на периодическую структуру подразумевает проведение нескольких численных экспериментов FDTD, которые мы называем итерациями. 
 +В течение каждой итерации происходит запись значений полей на периодических границах. 
 +Это дает ключ к решению проблемы получения полей из будущего: поля в будущем известны из прошлой итерации, поскольку они были уже записаны. 
 +При этом поля из прошлого могут браться из текущей итерации. 
 +Для управления итерационным процессорм используется метод Total Field/Scattered Field (TF/SF). 
 +Спустя некоторое число итераций граница TF/SF работает как применение периодических граничных условий с нужным сдвигом во времени. 
 + 
 +Наш метод описан в статье 
 +[[http://www.opticsinfobase.org/abstract.cfm?uri=ol-33-13-1491|http]]{{:valuev_-_iterative_technique_for_analysis_of_periodic_structures_at_oblique_incidence_in_the_finite-difference_time-domain_method.pdf|PDF}} 
 + 
 +{{tutorial_periodic:geometry_obl.png?350}}\\ 
 + 
 +Проиллюстрируем работу итерационного метода на примере расчета наклонного падения на фотонно-кристаллическую пластнику из упакованных в квадратную решетку сфер. 
 +Мы используем метод Total Field/Scattered Field (TF/SF) для генерации плоской волны внутри области Total Field, которая заштрихована на рисунке.  
 +TF/SF граничные условия на границе 3 задают наклонно падающую волну, которая известна заранее.  
 +Идея итерационного метода заключается в применении дополнительной TF/SF генерации на границах <tex>{\bf x}_1</tex> и <tex>{\bf x}_2</tex>, используя записанные значения полей из точек-образов 
 +<tex>{\bf x}_{1'}={\bf x}_1+{\bf a}</tex> и <tex>{\bf x}_{2'}={\bf x}_2-{\bf a}</tex>.  
 +Для отрицательного временного сдвига (поля из прошлого) можно брать поля из текущей итерации <tex>i</tex>, тогда как для положительного временного сдвига (поля из будещего) мы используем поля, записанные на предыдушей итерации: 
 + 
 +<tex> 
 +{\bf F}_i({\bf x}_2,t)={\bf F}_i({\bf x}_{2'},t-\Delta t), \\ 
 + 
 +{\bf F}_i({\bf x}_1,t)={\bf F}_{i-1}({\bf x}_{1'},t+\Delta t). 
 +</tex> 
 + 
 +Расстояние между границами 1 и 2 больше периода <tex>a</tex> на некоторую величину <tex>\Delta x</tex> равную нескольким сеточным шагам. Это нужно для разделения границ TF/SF 1, 2 и точек-образов.  
 + 
 +На первой итерации "падающая волна" на границах TF/SF 1, 2 есть нуль. 
 + 
 +Более подробно про итерационный метод можно прочитать здесь 
 +[[http://www.opticsinfobase.org/abstract.cfm?uri=ol-33-13-1491|http]]{{:valuev_-_iterative_technique_for_analysis_of_periodic_structures_at_oblique_incidence_in_the_finite-difference_time-domain_method.pdf|PDF}} 
 + 
 +Для включения итерационного метода нужно использовать функцию SetOblique. 
 +Первый параметр этой функции есть число итераций. 
 +Второй параметр (назовем его n) отвечает за работу детекторов: детекторы работают и анализируется на каждой n-ой и на последней итерации. 
 +Текстовые файлы детекторов с промежуточных итераций помечаются суффиксом  "_i" (i это номер итерации) после имени детектора. 
 + 
 +Ниже мы приводим код для моделирования наклонного падения в вакууме, рассматривается двумерный случай (размер вычислительного объема в направлении x есть нуль). 
 + 
 +<code cpp> 
 +#include "uiexp.h" 
 + 
 +int main(int argc,char **argv){ 
 + 
 +  emInit(argc,argv); 
 +  uiExperiment task; 
 + 
 +  task.SetInternalSpace(Vector_3(.5,0,0),Vector_3(.5,2,3)); 
 +  task.SetResolution(.05); 
 +  task.SetBC(BC_PER,BC_PER,BC_PML); 
 + 
 +  theta=45*M_PI/180.; // angle of incidence, radians. theta=0 corresponds to normal incidence 
 +  Vector_3 k=Vector_3(0,sin(theta),cos(theta)); // light propagation direction 
 +  Vector_3 E(1,0,0); // s polarization 
 +//  Vector_3 E(0,cos(theta),-sin(theta)); // p polarization 
 + 
 +  task.AddTFSFPlane(INF,INF,.5); 
 +  task.SetPlaneWave(k,E); 
 +  task.SetOblique(10,1); // comment it for normal incidence 
 + 
 +  task.BuildDimensions(); 
 + 
 +  Vector_3 p1,p2; 
 +  task.GetSpace(p1,p2); 
 +  Vector_3 mv1=p1+.05, mv2=p2-.05; // one mesh step from the calculated volume border 
 +  mv1[0]=mv2[0]=(p1[0]+p2[0])/2; 
 +  task.AddDetectorSet("movie",mv1,mv2,INT_INFTY,DET_TSTEP); 
 + 
 +  task.CalculateN(300); 
 + 
 +  task.Analyze(); 
 + 
 +  return 0; 
 +
 + 
 +</code> 
 + 
 +Заметим, что фактический размер вычислительного объема больше чем тот, который задан в аргументах функции SetInternalSpace. 
 +Этот размер включает в себя PML и пространство слева от границы 1 и справа от границы 2. 
 +Фактический размер вычислительного объема можно получить с помощью функции GetSize. 
 +Для иллюстрации сходимости метода, мы поставили детекторы "movie" по всему вычислительному объему. 
 +В течение численного эксперимента мы получим множество файлов moviexx_tyyyy.d, где xx это номер итерации, а yyyy это временной шаг. 
 +Эти файлы можно использовать для построения эволюции распространения поля на разных итерациях. 
 + 
 +  gnuplot> set pm3d interpolate 2,2 
 +  gnuplot> cbr=1.3 
 +  gnuplot> set palette defined (-cbr 'blue', 0 'white', cbr 'red') 
 +  gnuplot> set cbrange [-cbr:cbr] 
 +  gnuplot> set view 0,0 
 +  gnuplot> sp 'movie01_t0160.d' u 3:4:5 w pm3d 
 + 
 +{{tutorial_periodic:iter_obl.png?750}}\\ 
 + 
 +Как можно видеть из рисунка, в нашем случае распределение поля на 160 временном шаге сходится к решению (плоский импульс в вакууме) за первые 5 итераций. 
 +Для получения решения на более позднем временном шаге требуется больше итераций. 
 + 
 +=====Фотонно-кристаллическая пластинка===== 
 + 
 +Фотонно-кристаллическая пластинка это структура конечная вдоль одного направления и периодическая вдоль остальных направлений. 
 +Ниже мы моделируем наклонное падение на такую пластинку состояющую из металлических (<tex>\epsilon=4</tex>, <tex>\sigma=2</tex>) сфер упакованных в квадратную решетку.  
 +Период решетки <tex>a=1</tex>, радиус сфер <tex>R=0.375</tex>, сеточный шаг равен 0.025. 
 + 
 +<code cpp> 
 +#include "uiexp.h" 
 + 
 +int main(int argc,char **argv){ 
 + 
 +  emInit(argc,argv); 
 + 
 +  uiExperiment task; 
 +  task.SetInternalSpace(0,Vector_3(1,1,3)); 
 +  task.SetResolution(.025); 
 +  task.SetBC(BC_PER,BC_PER,BC_PML); 
 + 
 +  valtype theta=45*M_PI/180.; // angle of incidence, radians. theta=0 corresponds to normal incidence 
 +  Vector_3 k(0,sin(theta),cos(theta)); 
 +  Vector_3 E(1,0,0); // s polarization 
 +//  Vector_3 E(0,cos(theta),-sin(theta)); // p polarization 
 + 
 +  task.AddTFSFPlane(INF,INF,.5); 
 +  task.SetPlaneWave(k,E); 
 +  task.SetOblique(100,1); 
 + 
 +  task.AddObject(emMedium(4,2),new Sphere(.375,Vector_3(.5,.5,1.5))); 
 + 
 +  task.AddDetectorSet("h",Vector_3(.5,.5,0.25),Vector_3(.5,.5,2.75),iVector_3(1,1,2)); 
 +  task.AddRTASet("flux",2,0.25,2.75); 
 + 
 +  task.CalculateN(2000); 
 + 
 +  task.Analyze(); 
 + 
 +  return 0; 
 +
 +</code> 
 + 
 +Посмтроим геометрию в gnuplotе: 
 + 
 +  gnuplot> sp 'cell.pol' w l, 'med.pol' w l, 'pmlo.pol' w l, 'pmli.pol' w l, \ 
 +  gnuplot> 'tf.pol' w l, 'sg.vec' w vec, 'itf.pol' w l, 'unit.pol' w l 
 + 
 +{{tutorial_periodic:obl_gnu.png?350}}\\ 
 + 
 +Согласно картинке выше, 'itf.pol' соответствует границам 1 и 2, а 'unit.pol' соответствует периодической ячейке ширины a (она ограничена границей между 1 and 2' и границей между 1' и 2). 
 + 
 +Сравним результаты для спектра прохождения, полученные с помощью FDTD и полуаналитического метода layer Korringa-Kohn-Rostoker (LKKR). 
 + 
 +{{tutorial_periodic:obl_pc_trans.png?450}}\\
  
-======Прохождение через фотоннокристаллическую пластинку======+Как можно видеть из сравнения, результаты для больших частот сходятся быстрее. 
 +При этом 5-10 итераций достаточно для получения спектра прохождения с нормальной точностью. 
 +Для окончательной сходимости метода в нашем случае требуется порядка 100 итераций. 
 +Вообще, число итераций, необходимое для сходимости метода, пропорционально числу временных шагов в численном эксперименте (в нашем примере это 2000).
  
 ======Радиационное затухание источника====== ======Радиационное затухание источника======
  
  
 
/home/kintechlab/fdtd.kintechlab.com/docs/data/attic/ru/tutorial.1330590173.txt.gz · Последние изменения: 2012/03/01 12:22 — deinega     Наверх