Диссертация (1105278), страница 7
Текст из файла (страница 7)
Приведённая выше методика не предполагает учёта влияния симметрии навид компонент матрицы поляризуемости, все компоненты вычисляются независимо друг от друга.Заметим также, что описанный способ нахождения матрицы поляризуемости (1.13) пригоден не только, когда значения комплексных амплитуд дипольных моментов получены с помощью аналитических или численных расчётов,но и в тех случаях, когда они измерены экспериментальными методами.2.3.Численное решение задачи рассеяния методом конечных элементов в COMSOL MultiphysicsКак было отмечено выше, в настоящей работе для решения задачи рас-сеяния использовался RF модуль пакета программ COMSOL Multiphysics1 , воснове которого лежит метод конечных элементов [103, 104].Рассмотрим основные этапы решения задачи рассеяния методом конечных элементов в COMSOL Multiphysics на примере П-образного резонатора,представленного на рис.
2.2. Форма резонатора задавалась следующими геометрическими параметрами: длинной стороны l, высотой резонатора h, ширинойзазора g и глубиной зазора u. В этой главе ориентация П-образного резонаторавдоль осей декартовых координат принята следующая: плоскость резонатораперпендикулярна оси y, ось z проходит вдоль ножек резонатора, а ось x —перпендикулярно им.1v.3.5a, build 3.5.0.60840Рис. 2.2: Геометрические параметры отдельного резонатора и его ориентация в пространстве.Частицу окружала сферическая область радиусом 2,5l, электромагнитныехарактеристики которой соответствуют вакууму.Для того, чтобы избежать паразитных отражений вакуумная сферическойобласть окружена областью представлявшей идеальную согласованную нагрузку (в англоязычной литературе используется термин PML, Perfectly matchedlayer), эффективно поглощающую рассеянное излучение и прошедшую волну.Толщина слоя PML составляла 0,5l.В области расчёта задаются материальные параметры веществ и граничные условия, сшивающие поля, а также амплитуда, направление и длина падающей плоской волны.В области расчёта находилось решение волнового уравнение для рассеянного электрического поля Esc :1σ2~ sc − k0 ε − i~ sc = 0.∇×EE∇×µωε0(2.38)Для этого вся область разбивалась на сетку, состоящую из множества конечных элементов в форме тетраэдров и призм, как показано на рис.
2.3.После разбиения области на конечные элементы COMSOL Multiphysicsпроизводит расчёт в автоматическом режиме без участия пользователя. Дляполноты изложения рассмотрим основные шаги метода конечных элементовна этапе собственно решения. В каждом конечном элементе неизвестное решение волнового уравнения (2.38) для электрического поля представлялось ввиде полиномиальной функции второй степени с неизвестными коэффициентами. Функции не равны нулю только внутри своих конечных элементов. Затем,методом Ритца или Галёркина, производился переход к системе (разреженных)линейных уравнений, решение которой давало искомое распределение поля в41Рис. 2.3: Область расчёта, разбитая на конечные элементы.
Показаны конечные элементы награницах областей.области расчёта. Подробное и детальное рассмотрение метода конечных элементов и его применение к задачам электродинамики можно найти в монографиях [103, 104].COMSOL Multiphysics предоставляет на выбор довольно большой наборразличных решателей (в англоязычной литературе широко используется термин solver) систем разреженных линейных уравнений. Довольно полное рассмотрение возможностей и сравнение производительности различных решателей приведено в работах [105,106]. На основании анализа этих работ и по результатам собственных тестовых расчётов, было определено, что наиболее эффективный для рассматриваемых проблем решатель — PARDISO. Он характеризуется довольно быстрой загрузкой оперативной памяти при старте моделирования и эффективным использованием многоядерных процессоров с поддержкойнескольких потоков.Первоначально для проведения расчётов использовалась ЭВМ оснащеннаяпроцессором Intel Core i7 950 (4 ядра, 8 потоков) и 12 Гб оперативной памяти.Среднее время расчёта распредления поля в представленной на рис.
2.3 областидля одной точки по частоте составляло от десятков секунд до десятков минут инелинейно зависело от количества элементов на которые была разбита областьрешения. Для модели состоящей из примерно 22 тысяч элементов среднее время42решения 107 секунд. Впоследствии расчёты также производились на машине спроцессорам Intel Core i7 3770k и 16 Гб оперативной памяти, что сократиловремя вычисления модели с 22 тысячами элементов до 70 секунд для однойточки по частоте.Возможность COMSOL Multiphysics работать в связке с MATLAB позволила написать ряд программ параметризующих и автоматизирующих процессмоделирования. Постобработка результатов расчёта в основном осуществляласьв MATLAB.После решения численной задачи по известному распределению полей всоответствии с формулами (2.22), (2.23) и (2.28), (2.29) рассчитываются компоненты индуцированных дипольных моментов.
Для полноты изложения приведём выражения, использующиеся при их вычислении в COMSOL Multiphysics,в явном виде.Объёмные компоненты электрического дипольного момента получаютсяпутём интегрирования по объёму частицы следующих выражений:epsilon0_rfw * (x * (d(Ex,x) + d(Ey,y) + d(Ez,z)))epsilon0_rfw * (y * (d(Ex,x) + d(Ey,y) + d(Ez,z)))epsilon0_rfw * (z * (d(Ex,x) + d(Ey,y) + d(Ez,z)))Поверхностные компоненты электрического дипольного момента получаются путём интегрирования следующих выражений на границах частицы:epsilon0_rfw * (x * (dnx*(up(Ex) - down(Ex)) ...+ dny*(up(Ey) - down(Ey)) ...+ dnz*(up(Ez) - down(Ez))))epsilon0_rfw * (y * (dnx*(up(Ex) - down(Ex)) ...+ dny*(up(Ey) - down(Ey)) ...+ dnz*(up(Ez) - down(Ez))))epsilon0_rfw * (z * (dnx*(up(Ex) - down(Ex)) ...+ dny*(up(Ey) - down(Ey)) ...+ dnz*(up(Ez) - down(Ez))))Объёмные компоненты магнитного дипольного момента получаются путёминтегрирования по объёму частицы следующих выражений:43i*pi*nu_rfw * ...(epsilon0_rfw*epsilonr_rfw - i*sigma_rfw/(2*pi*nu_rfw)) * ...(y*Ez - z*Ey)i*pi*nu_rfw * ...(epsilon0_rfw*epsilonr_rfw - i*sigma_rfw/(2*pi*nu_rfw)) * ...(z*Ex - x*Ez)i*pi*nu_rfw * ...(epsilon0_rfw*epsilonr_rfw - i*sigma_rfw/(2*pi*nu_rfw)) * ...(x*Ey - y*Ex)Поверхностные компоненты магнитного дипольного момента получаютсяпутём интегрирования следующих выражений на границах частицы:1/2 * (y*Jsz_rfw - z*Jsy_rfw)1/2 * (z*Jsx_rfw - x*Jsz_rfw)1/2 * (x*Jsy_rfw - y*Jsx_rfw)2.4.Дипольные поляризуемости сферических частицДля проверки работоспособности и правильности задания моделей внача-ле были рассчитаны поляризуемости сферических частиц.
С этой целью, кроме COMSOL Multiphysics, также были использованы пакеты программ AnsoftHFSS и CST Microwave Studio. Компоненты электрической и магнитной поляризуемости в дипольном приближении были определены с помощью формул(2.22), (2.23) и (2.28), (2.29). В силу симметрии, тензоры поляризуемости дляeeeшара содержат только диагональные компоненты αe = αxx= αyy= αzzиmmmαm = αxx= αyy= αzz.Эти параметры могут быть также получены исходя из аналитических выражений теории Ми [9, 107]. Сравнение результатов численного моделированияи расчётов по точным аналитическим формулам даст прямой и самый точныйспособ оценки величины погрешности, которая вносится при моделировании.Как было показано в работе [91, 108] диагональные компоненты определяются следующими соотношениями (в системе СИ):44αe = 6πia1 /k23 ,αm = 6πib1 /k23 ,(2.39)где a1 и b1 — первые парциальные амплитуды разложения рассеянного излучения в ряд по сферическим гармоникам.
Общие выражения для парциальныхамплитуд теории Ми записываются в виде [90]:mψn (mx)ψn0 (x) − ψn (x)ψn0 (mx)an =,mψn (mx)ξn0 (x) − ξn (x)ψn0 (mx)ψn (mx)ψn0 (x) − mψn (x)ψn0 (mx)bn =,ψn (mx)ξn0 (x) − mξn (x)ψn0 (mx)(2.40)где ψn и ξn — функции Риккати-Бесселя 1-го и 2-го рода (функции Бесселяот полуцелого индекса), m = ñ1 /ñ2 — комплексный относительный показательпреломления.Расчёт поляризуемости по формулам теории Ми был выполнен с помощью программы, разработанной в диссертации А.В.
Журавлева [107]. Сравнение результатов по теории Ми и численному моделированию был произведёндля случая рассеяния плоской электромагнитной волны на металлическом шаре диаметром 60 нм. В качестве материала частиц были выбраны алюминийи золото, как материалы с наиболее высокими значениями плазменной частоты.
Значения комплексной диэлектрической проницаемости, необходимые длязадания модели, взяты из экспериментальной работы [109]. Частотные зависимости реальной части диэлектрической проницаемости, проводимости и толщины согласно данным из [109] для золота и алюминия представлены на рис. 2.4.Дополнительно приведены спектры для золота согласно работе [110], котораятакже широко используется в литературе.Результаты сравнения расчётов по аналитическим формулам теории Мии методами численного моделирования представлены на рис. 2.5.
Расчёты показали согласие зависимостей поляризуемостей, вычисленных разными методами. Относительная погрешность моделирования не превышала 5%. Следуетотметить, что при увеличении отношения радиуса частицы к длине волны еёизлучение в данном спектральном диапазоне перестаёт быть дипольным и использование формул (2.22), (2.23) и (2.28), (2.29) даёт погрешности порядка 20%Границы областей, где излучение частицы можно считать дипольным, можно45Рис. 2.4: Частотные зависимости реальной (а) и мнимой (б) частей диэлектрической проницаемости и толщины скин слоя (в) золота и алюминия (эксперимент). Сплошные линии —алюминий [109], пунктир — золото [109], точечные линии — золото [110] — точечная линия.определить с помощью графиков отношения мощности дипольного излученияк полной мощности излучения, рассчитанных в рамках теории Ми (см.
рис.2.6) [107]. Белые линии D(λ) на рис. 2.6 соответствуют уровню в 95%. Для образцов меньших размеров (диаметр d < D(λ)) дипольное приближение можносчитать примениным. На рис. 2.5 область, где справедливо дипольное приближение соответствует длине волны большей 500 нм.На основании результатов моделирования можно сделать вывод о пригодности COMSOL Multiphysics для вычисления поляризуемостей резонаторов.Другие программные пакеты обеспечивают сопоставимые погрешности и скорости расчёта. COMSOL Multiphysics был выбран из-за более логично построенного и интуитивного интерфейса (на взгляд автора), удобных инструментовпост-процессинга и глубокой интеграции со средой MATLAB. Последнее позволило, путём написания ряда программ на MATLAB, параметризовать и автоматизировать процесс расчёта больших и ресурсоёмких моделей.Аналогичный вывод можно сделать, рассматривая работу [111], в которой использовалась схема моделирования довольно близкая к использованнойв настоящей работе и получено хорошее соответствие сечений рассеяния сферической частицы, полученных моделированием в COMSOL Multiphysics, с сечениями рассеяния рассчитанными по теории Ми.46Рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
