Различия

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

ru:tutorial [2012/07/22 22:02]
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>
  
 =====Детекторы===== =====Детекторы=====
Строка 576: Строка 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).
 Пластинка сделана из дисперсной среды, описываемой двумя членами Лоренца. Пластинка сделана из дисперсной среды, описываемой двумя членами Лоренца.
  
Строка 594: Строка 569:
  
 Поглощение может быть получено как 1 - R - T. Поглощение может быть получено как 1 - R - T.
 +
 +Расстояние между телами, границей TF/SF, PML и детекторами должно быть не менее 2 сеточных шагов. В нашем примере мы придерживались большего расстояния (5-20 сеточных шагов), поскольку иногда это может улучшить результат.
  
 =====Кремниевая подложка===== =====Кремниевая подложка=====
Строка 599: Строка 576:
 Ниже мы рассчитываем отражение от полубесконечной кремниевой подложки в случае нормального падения. Мы используем уже имеющийся код из предыдущей подсекции (задание граничных условий, области Total Field и т. д.) Ниже мы рассчитываем отражение от полубесконечной кремниевой подложки в случае нормального падения. Мы используем уже имеющийся код из предыдущей подсекции (задание граничных условий, области Total Field и т. д.)
  
-Для задания полубесконечной подложки мы можем воспользоваться функцией GetPolyhedronPlane. В отличие от пластинки, у подложки нет второй замыкающей поверхности, и она бесконечно продолжается вдоль оси z. В модели FDTD она погружается в PML, что моделирует ее бесконечность.+Для задания полубесконечной подложки мы можем воспользоваться функцией GetHalfSpace. В отличие от пластинки, у подложки нет второй замыкающей поверхности, и она бесконечно продолжается вдоль оси z. В модели FDTD она погружается в PML, что моделирует ее бесконечность.
  
-Материалом подложки будет кремний. Для получения соответствующего объекта emMedium можно воспользоваться функцией getSi, у которой имеется аргумент, по умолчанию равный единице. Этот аргумент задает, чему равна единица длины FDTD в микронах. Это необходимо для правильной калибровки коэффициентов, участвующих в расчете диэлектрической проницаемости (см. [[tutorial?&#media|see how to specify media in EMTL]]). В нашем случае единица FDTD будет равна 0.1 микрона.+Материалом подложки будет кремний. Для получения соответствующего объекта emMedium можно воспользоваться функцией getSi, у которой имеется аргумент, по умолчанию равный единице. Этот аргумент задает, чему равна единица длины FDTD в микронах. Это необходимо для правильной калибровки коэффициентов, участвующих в расчете диэлектрической проницаемости (см. [[tutorial?&#sredy|как задавать среды в EMTL]]). В нашем случае единица FDTD будет равна 0.1 микрона.
  
-  task.AddObject(getSi(.1),GetPolyhedronPlane(Vector_3(0,0,1), Vector_3(0,0,1)));+  task.AddObject(getSi(.1),GetHalfSpace(Vector_3(0,0,1), Vector_3(0,0,1)));
      
 Параметры, используемые для апроксимации диэлектрической проницаемости кремния функцией getSi, были получены с помощью программы на MatLab, которую вы можете найти в разделе [[fitting]]. Параметры, используемые для апроксимации диэлектрической проницаемости кремния функцией getSi, были получены с помощью программы на MatLab, которую вы можете найти в разделе [[fitting]].
Строка 619: Строка 596:
  
 Как видно, кремний достаточно сильно отражает в оптическом диапазоне. Для уменьшения отражения могут использоваться [[ar|антиотражающие текстурированные покрытия]]. Как видно, кремний достаточно сильно отражает в оптическом диапазоне. Для уменьшения отражения могут использоваться [[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.1342980145.txt.gz · Последние изменения: 2012/07/22 22:02 — deinega     Наверх