На правах рукописи

Дейнега Алексей Вадимович

Численное моделирование и компьютерный дизайн

оптических свойств

наноструктурированных материалов

Специальность 01.04.05 Оптика

Автореферат

диссертации на соискание ученой степени

кандидата физико-математических наук

Москва 2010

Работа выполнена в Институте водородной энергетики и плазменных технологий Российского Научного Центра Курчатовский Институт.

Научный руководитель: кандидат физико-математических наук, доцент Потапкин Борис Васильевич

Официальные оппоненты: доктор физико-математических наук, доцент Ёлкин Николай Николаевич (Троицкий институт инновационных и термоядерных исследований, г. Троицк) доктор физико-математических наук, профессор Козарь Анатолий Викторович (Физический факультет МГУ имени М. В. Ломоносова, г. Москва)

Ведущая организация: Институт спектроскопии РАН (г. Троицк)

Защита состоится 19 мая 2010 г. в 15 часов на заседании совета по защите докторских и кандидатских диссертаций Д 501.001.45 при Московском Государственном Университете имени М. В. Ломоносова по адресу: Россия, 119991, Москва, Ленинские горы, дом 1, строение 5 ( корпус НИИ ядерной физики МГУ), ауд. 2-15.

С диссертацией можно ознакомиться в научной библиотеке НИИ ядерной физики имени Д. В. Скобельцына МГУ имени М. В. Ломоносова.

Автореферат разослан 15 апреля 2010 г.

Ученый секретарь совета по защите докторских и кандидатских диссертаций Д 501.001.45 при МГУ имени М. В. Ломоносова кандидат физико-математических наук Вохник О. М.

1

Общая характеристика работы

Диссертация посвящена численному моделированию оптических свойств наноструктурированных материалов, а именно, в ней исследуются следующие вопросы:

Возможность применения металлических фотонных кристаллов в качестве новых высокоэффективных источников света.

Оптимизация размеров и формы антиотражающих нанотекстурированных покрытий.

Основным используемым численным методом является метод решения уравнений Максвелла в конечных разностях (Finite-Dierence Time-Domain, FDTD). Для расчетов применяется специально написанная параллельная программа, включающая в себя ряд новых оригинальных численных методов, которые также описываются в данной диссертации.

Актуальность работы. Наноструктурированные материалы представляют собой новый тип материалов, характеризующийся малым размером (характерный размер порядка нескольких сотен нанометров и менее) и сложной организацией составляющих его элементов. Наноматериалы обладают уникальными физическими свойствами, что позволяет им находить применение во многих промышленных областях: в вычислительной технике, энергетике, медицине и т. д. [1] К наноматериалам относятся фотонные кристаллы [2], представляющие собой структуры с периодически меняющейся в пространстве диэлектрической проницаемостью. Оптические свойства фотонных кристаллов сильно отличаются от оптических свойств сплошных сред: электромагнитные волны в фотонном кристалле имеют зонный спектр и координатную зависимость, аналогичную блоховским волнам электронов в обычном кристалле. Перераспределение плотности фотонных состояний в фотонном кристалле приводит к изменению спектра их теплового излучения [3]. Внутри запрещенной зоны плотность фотонных состояний равна нулю, и излучение электромагнитных волн подавлено. Вне запрещенной зоны интенсивность этого излучения непосредственно связана с плотностью фотонных состояний и может сильно варьироваться на разных частотах. На основании этого в литературе было высказано предположение о возможности использования металлических фотонных кристаллов в качестве высокоэффективных источников видимого света [3], излучение которых было бы подавлено в инфракрасном диапазоне, что позволило бы добиться существенно большей световой отдачи по сравнению с обычными лампами накаливания. Однако до сих пор не было проведено какого-либо детального исследования практического осуществления такой возможности. В частности, не было изучено влияние неизбежных при изготовлении фотонных кристаллов дефектов, а также роли внешней матрицы, необходимой для фиксации составляющих фотонный кристалл элементов. В связи с этим становится актуальным исследование излучательных характеристик металлических фотонных кристаллов при высоких температурах и возможности создания на их основе высокоэффективных источников света.

Еще одним примером практического использования наноматериалов являются текстурированные антиотражающие покрытия. Интерес к ним особенно возрос в последнее время, а именно, появилась масса новых работ по успешному изготовлению антиотражающих нанотекстурированных покрытий (см., напр., [4]), используемых, в частности, в солнечной энергетике. Параллельно активно ведется теоретическое исследование оптических свойств текстурированных покрытий (см., напр., [5]). Отметим, что в имеющихся работах, как правило, либо отдельно рассматриваются длинно- или коротковолновый предельные случаи, либо видимый диапазон для размеров составляющих текстуру элементов порядка нескольких сотен нанометров. В это же время является актуальным единое исследование текстурированных поверхностей во всем диапазоне длин волн. Такое исследование позволило бы установить, каков должен быть размер составляющих текстуру элементов, чтобы для заданного диапазона длин волн отражение было бы минимальным.

В этой работе основным используемым численным методом является метод решения уравнений Максвелла в конечных разностях (англ. - FiniteDierence Time-Domain method, FDTD) [6]. Выбор этого метода вызван присущей ему высокой параллельной эффективностью, что позволяет применять его для расчета больших задач на кластерных вычислительных системах.

Как и всякий другой разностный метод, FDTD сталкивается с проблемой неточного отображения на прямоугольную вычислительную сетку тел, обладающих произвольной формой поверхности. Решение этой проблемы особенно значимо для моделирования исследуемых в данной работе фотонных кристаллов и текстурированных поверхностей. Одним из наиболее удачных способов ее решения является метод подсеточного сглаживания [7]. Однако этот метод в его первоначальной постановке может применяться только для диэлектриков, в связи с чем актуальна его модификация, которая позволила бы расширить его на дисперсные среды.

Исследование оптических свойств наноматериалов включает в себя получение спектров прохождения и отражения от периодических структур при падении на них плоской волны. Их расчет FDTD для случая нормального падения предполагает применение периодических граничных условий. В случае наклонного падения периодические граничные условия содержат временной сдвиг, что затрудняет их применение. В литературе было предложено несколько методов для решения этой проблемы, однако у каждого из них есть свои недостатки: нестабильность для углов близких к, увеличение объема используемой памяти и т. д. В связи с этим актуальна разработка альтернативного метода, который был бы лишен недостатков предыдущих. Такой метод позволил бы эффективно исследовать излучательные характеристики фотонных кристаллов, а также антиотражающие свойства текстурированных покрытий при любом угле падения.

Цель работы состоит в 1) исследовании излучательных характеристик металлических фотонных кристаллов, а также возможности их практического применения в качестве источников света c высокой селективностью излучения в видимом диапазоне; 2) исследовании антиотражающих свойств текстурированных покрытий во всем диапазоне размеров составляющих текстуру рассеивателей, а также установлении оптимального размера для заданного диапазона длин волн; 3) дальнейшем развитии метода FDTD, включающем, в частности, разработку метода подсеточного сглаживания для дисперсных сред и метода расчета наклонного падения плоской волны на периодическую структуру.

Научная новизна работы.

1. Впервые исследованы излучательные характеристики вольфрамовых фотонных кристаллов типа прямого опала, погруженных в керамическую матрицу. Найдены геометрические параметры фотонного кристалла (период решетки порядка нескольких сотен нанометров, фактор заполнения по вольфраму порядка нескольких процентов), при которых в случае отсутствия поглощения во внешней матрице КПД источника света на основе фотонного кристалла оказывается выше в несколько раз по сравнению с обычными лампами накаливания. Показано, что у рассмотренных составляющих матрицу термоустойчивых материалов (оксид гафния) сильное поглощение при высоких температурах приводит к существенному снижению эффективности светового источника на основе фотонного кристалла, вплоть до значения КПД лампы накаливания. Сделан вывод о том, что для создания высокоэффективного источника света на основе фотонного кристалла необходим поиск материалов со слабым поглощением при высоких температурах.

2. Впервые исследовано влияние дефектов в фотонном кристалле (неупорядоченности по расположению и разброса по размерам составляющих фотонный кристалл элементов) на его излучательную способность. Показано, что их наличие приводит к уменьшению величины излучения в видимом диапазоне при сохранении подавления излучения в инфракрасном диапазоне.

Получено количественное согласие результатов расчетов с результатами специально поставленного эксперимента.

3. Впервые проведено систематическое сравнительное исследование антиотражающих свойств текстурированных покрытий во всем диапазоне размеров составляющих текстуру рассеивателей по отношению к длине волны.

Установлено, что ключевым фактором, влияющим на оптимальный размер составляющих текстуру рассеивателей, является характер замощения подложки их основаниями. В случае полного замощения наименьшее значение отражения достигается при макроскопических размерах рассеивателей, а в случае неполного - оптимальный размер имеет порядок длины волны.

4. Впервые установлены асимптотические зависимости величины отражения от геометрических параметров текстурированной поверхности для длиннои коротковолнового пределов. В случае полного замощения подложки основаниями составляющих текстуру рассеивателей в длинноволновом пределе отражение уменьшается с увеличением высоты рассеивателей степенным образом, а в коротковолновом - экспоненциальным; в случае неполного замощения в длинноволновом и коротковолновом пределах величина отражения выходит на постоянное значение.

5. Предложен новый подход программной реализации метода решения уравнений Максвелла в конечных разностях (FDTD), основанный на интегральном представлении уравнений Максвелла. В рамках этого подхода написана параллельная программа на языке программирования С++. В ходе ее написания реализован ряд оригинальных решений: оптимизация использования памяти с помощью упаковки используемых в разностных уравнениях коэффициентов; увеличение производительности посредством разбиения алгоритма обновления сеток на основной цикл, поддающийся векторизации, и специализированные поправки к нему; увеличение параллельной эффективности путем балансировки доменов. Специально проведенные тесты продемонстрировали линейный характер масштабируемости программы вплоть до тысячи вычислительных ядер. Программа находится в открытом доступе в интернете (http://fdtd.kintechlab.com).

6. В рамках метода FDTD впервые предложен и реализован алгоритм подсеточного сглаживания для дисперсных сред, позволяющий существенно уменьшить величину численной ошибки, связанной с неточным отображением тел, обладающих произвольной формой поверхности, на прямоугольную вычислительную сетку. Его применение приводит к увеличению порядка точности метода FDTD по сравнению с обычным лестничным приближением.

7. Впервые предложен и реализован итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру в численном расчете FDTD. С его помощью можно получать результаты во всем диапазоне углов от до, а также наблюдать реальное поведение поля во времени. Помимо этого, он лишен таких недостатков предложенных ранее методов, как неустойчивость при отвесных углах падения и увеличение требуемого размера сетки. Особенностью метода является то, что получение решения для одного угла требует нескольких отдельных итераций-экспериментов.

Положения, выносимые на защиту.

1. Геометрические параметры вольфрамовых фотонных кристаллов типа прямого опала, при которых КПД источника света на основе фотонных кристаллов выше в несколько раз по сравнению с обычными лампами накаливания.

2. У рассмотренных термоустойчивых материалов (оксид гафния), используемых в матрице, необходимой для фиксации элементов фотонного кристалла, сильное поглощение при высоких температурах приводит к существенному снижению КПД светового источника на основе фотонного кристалла, вплоть до значения КПД лампы накаливания.

3. Асимптотические зависимости величины отражения от геометрических параметров текстурированной поверхности для длинно- и коротковолнового пределов.

4. Ключевым фактором, влияющим на оптимальный размер составляющих текстуру рассеивателей, является характер замощения подложки их основаниями: в случае полного замощения наименьшее значение отражения достигается при макроскопических размерах рассеивателей, а в случае неполного - оптимальный размер рассеивателей имеет порядок длины волны.

5. Параллельная программа FDTD, которая может использоваться для предсказательного моделирования в разнообразных оптических приложениях.

6. Новые численные методы в рамках FDTD: метод подсеточного сглаживания для дисперсных сред и итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру.

Практическая ценность работы.

1. Найдены геометрические параметры вольфрамовых фотонных кристаллов типа прямого опала, при которых КПД источника света на основе фотонного кристалла оказывается выше в несколько раз по сравнению с обычными лампами накаливания.

2. Установлено, что для создания высокоэффективных источников света на основе фотонного кристалла необходим поиск материалов со слабым поглощением при высоких температурах.

3. Найдена оптимальная геометрия текстурированной поверхности, при которой достигается наименьшая величина отражения для заданного диапазона длин волн.

4. Показано, что для достижения малых значений отражения в видимом диапазоне наиболее эффективным решением являются нанотекстурированные покрытия.

5. Разработана параллельная программа FDTD, которая может использоваться для предсказательного моделирования в самых разнообразных оптических приложениях: при моделировании оптических свойств фотонных кристаллов, при расчете антиотражающих свойств текстурированных покрытий и т. д. Программа находится в открытом доступе в интернете.

6. Предложены и реализованы численные методы, позволяющие существенно увеличить быстродействие и точность расчетов FDTD.

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

Автором получены основные результаты для излучательных характеристик вольфрамовых фотонных кристаллов типа прямого опала.

Автор внес решающий вклад при написании параллельной программы, реализующей метод FDTD, а также при реализации итерационного метода для расчета наклонного падения плоской волны на периодическую структуру.

Автором был самостоятельно предложен и реализован метод подсеточного сглаживания для дисперсных сред.

Достоверность результатов. Достоверность расчетов обеспечивается системой многократных проверок имеющегося кода, а также путем сравнения полученных с помощью него результатов с результатами, полученными другими методами и с имеющимися в литературе экспериментальными данными.

Апробация работы. Основные результаты диссертационной работы были доложены на научно-технических конференциях The International Conference on Computational Science and Its Applications (Малайзия, Kuala Lumpur, 2007), Conference on Computational Physics (Бразилия, Ouro Preto, 2008), Conference on Computational Physics (Тайвань, Kaohsiung, 2009), Многомасштабное моделирование процессов и структур в нанотехнологиях (Москва, 2009), VII Курчатовская молодежная научная школа (Москва, 2009). 52-ая научная конференция МФТИ (Москва, 2009).

Реализованная программа FDTD стала призером на конкурсе Максимальная масштабируемость, проведенном в рамках РОСНАНОФОРУМАкомпаниями Intel и "Роснано".

Публикации. По материалам диссертации опубликовано 5 работ в реферируемых научных изданиях, 1 работа в нереферируемом научном журнале и тезисы российских и международных конференций (список в конце автореферата).

Структура и объем диссертации. Диссертация состоит из введения, четырех глав и заключения, изложена на 157 страницах, включая 74 рисунков и 125 наименований цитируемой литературы.

2 Содержание работы Во введении обоснована актуальность проводимых в работе исследований, сформулированы основные задачи диссертационной работы, оценена их научная новизна, изложена структура диссертации.

В первой главе обсуждается метод решения уравнений Максвелла в конечных разностях (Finite-Dierence Time-Domain, FDTD). Эта глава содержит в себе описание метода FDTD и изложение новых предложенных численных методов: метод подсеточного сглаживания для дисперсных сред и итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру.

Метод подсеточного сглаживания предназначен для уменьшения численной ошибки, связанной с неточным отображением тел с произвольной формой поверхности на вычислительную сетку. Он основан на введении тензора обратной эффективной диэлектрической проницаемости вблизи границы между телами [7]:

где - это матрица соответствующая проектору на вектор нормали к границе между двумя средами, а есть усреднение по контрольному объему, окружающему выбранный узел сетки.

До сих пор метод подсеточного сглаживания применялся только для моделирования диэлектриков. Мы развили его также и для случая дисперсивных сред. Наша идея заключается в разбиении электрического поля на вспомогательные компоненты таким образом, что получающиеся уравнения для каждой компоненты могут быть решены во временном представлении традиционным для FDTD путем. Например, если зависимость диэлектрической проницаемости от частоты представляется в виде суммы произвольного числа членов Друде или Лоренца, эти уравнения могут быть решены методом вспомогательных дифференциальных уравнений [6].

Для проверки работы изложенного метода мы сравнили полученные с помощью него результаты для рассеяния плоской волны на сфере с аналитическим решением для этой задачи. Мы проводили тестирование следующих efficiency factor for scattering Рис. 1: Слева: Зависимость фактора эффективности рассеяния от радиуса для свинцовой сферы,.

Справа: Угловая зависимость элемента матрицы рассеяния для свинцовой сферы радиуса nm,.

методов: 1) лестничное приближение, в котором значение диэлектрической проницаемости берется в узле сетки; 2) усреднение диэлектрической проницаемости по контуру или по объему; 3) усреднение обратной диэлектрической проницаемости; 4) VP-EP метод [7], который был нами реализован для дисперсивных сред; 5) собственно применение полного выражения (1).

Рассмотрим результаты расчетов для рассеяния плоской волны на свинцовой сфере (радиус варьируется), рассчитанное на сетке с шагом nm (Рис. 1). Отметим, что на рассматриваемой длине волны действительная часть диэлектрической проницаемости свинца отрицательна:

. Поскольку решение уравнений Максвелла существенно различается для областей с положительной и отрицательной диэлектрической проницаемостью, методы, не учитывающие направление границы между телами и основанные на частичной аппроксимации выражения (1), показывают существенное снижение точности по сравнению новым методом, в котором примененяется тензор обратной эффективной диэлектрической проницаемости.

Далее в первой главе мы предлагаем новый итерационный метод для моделирования наклонного падения плоской волны на периодическую планарную структуру. Наш метод предполагает проведение нескольких ( - ) отдельных численных экспериментов FDTD, которые мы в дальнейшем будем называть итерациями. Идея метода заключается в итерационном приближении к искомому решению, для которого поле в граничных точках периодической ячейки и в любой момент времени удовлетворяет условию:

Рис. 2: Геометрия численного эксперимента. 1, 2 – генерирующие (TF/SF) границы, для которых в качестве источника берется сигнал, переносимый с точек образов, ; 3 – TF/SF граница для генерации падающей плоской волны; 4 – детекторы для записи прошедшего сигнала.

где есть периоды структуры по двум направлениям, - величина угла падения в плоскости, образованной соответствующим направлением трансляции двумерной решетки и нормалью к ней. В зависимости от положения рассматриваемой границы ячейки применяемый временной сдвиг может быть как положительным, так и отрицательным, а, следовательно, значения полей в противоположных приграничных точках должны браться либо из прошлого, либо из будущего.

Геометрия численного эксперимента, соответствующего предлагаемому нами методу, представлена на Рис. 2 (в целях простоты изложения мы рассматриваем двумерный случай с, и временным сдвигом ). Мы используем метод общего и рассеянного поля (Total Field / Scattered Field Method, TF/SF) [6] для генерации плоской волны в область общего поля (эта область заштрихована на Рис. 2). Зависящие от времени граничные условия TF/SF на границе 3 соответствуют падающей наклонно плоской волне и известны аналитически. Основная идея нового метода заключается в применении дополнительных граничных условий на границах 1 и 2, которые работают подобно границе TF/SF, где в качестве генерируемой волны используются значения полей в точках-образах на границах Для отрицательного временного сдвига в прошлое могут быть использованы значения полей с текущей -ой итерации, а для положительного сдвига в будущее мы используем значения полей, записанные на предыдущей -ой итерации:

Transmittance Рис. 3: Слева: Зависимость коэффициента прохождения однослойного фотонного кристалла, состоящего из упакованных в квадратную решетку металлических шариков (, Справа: Зависимость коэффициента прохождения от угла падения, частота.

торое расстояние, равное нескольким сеточным шагам, для того чтобы разделить точки-образы и дополнительную границу TF/SF. Для записи значений и на границе 1 используется дополнительный буфер памяти размером, где - продолжительность одной итерации, - шаг по времени, а - количество сеточных узлов, приходящихся точках-образах полагаются равными нулю. На первой итерации в качестве граничного условия на границах 1 и 2 берется аналитическое значение для падающей плоской волны или нулевой сигнал.

Далее в первой главе приводится доказательство сходимости предложенного метода. Идея доказательства заключается в представлении генерируемого на -ой итерации на границах 1 и 2 сигнала в виде суммы, где удовлетворяет условию (2), а есть сигнал-ошибка. По причине линейности уравнений Максвелла и операции временного сдвига по отношению к полям, эволюция ошибки при увеличении может рассматриваться отдельно от решения. Нами доказано и численно проверено, что сигнал-ошибка затухает при увеличении числа итераций.

Метод был протестирован при моделировании однослойного фотонного кристалла, состоящего из упакованных в квадратную решетку металлических шариков (, ). Период квадратной решетки полагался равным единице, радиус сферы, а шаг по сетке 0.05. На Рис. представлено сравнение результатов для спектра прохождения, полученных с помощью описанного метода и с помощью метода LKKR [8]. Результаты FDTD находятся в хорошем согласии с результатами LKKR уже после пяти итераций.

Новый метод лишен таких недостатков предложенных раннее в литературе методов [6], как неустойчивость при отвесных углах падения, увеличение требуемого размера сетки, переход в частотное представление.

Во второй главе описывается использовавшийся нами подход для написания параллельного кода FDTD. В самом широком смысле метод FDTD включает в себя множество самых разнообразных возможностей: применение различных типов сеток, использование методов постпроцессорной обработки результатов и т. д. Для реализации всех этих возможностей требуется написание больших программ, которые обладали бы высокой производительностью, были бы легко расширяемы и обладали доступным, читаемым кодом.

К сожалению, распространенные некоммерческие программы FDTD лишены этих свойств и не предназначены для широкого круга задач. Это явилось причиной разработки собственной модели FDTD, которую мы будем в дальнейшем называть EMTL (Electromagnetic Template Library). EMTL реализована в виде библиотеки шаблонов C++, которые могут быть использованы в качестве строительных компонент при создании конкретных численных приложений. Для реализации возможности проведения параллельных расчетов используется интерфейс MPI (Message Passing Interface).

EMTL можно представить как совокупность взаимодействующих друг с другом компонент. Базовые компоненты могут быть задействованы в качестве параметров шаблонов при построении компонент более высокого уровня.

Такими базовыми компонентами являются тела, итераторы, сетки и т. д.

Ключевым строительным блоком EMTL является вычислительная сетка.

Сетка организует размещение данных в памяти и определяет способ интерполяции значений поля по сеточным узлам. В рамках этих требований могут быть реализованы сетки самого разного типа, включая не ортогональные и неоднородные. Проведение численного расчета с использованием сетки и внутреннее устройство самой сетки рассматриваются в EMTL как две независимые задачи. Это позволяет реализовать многие процедуры EMTL в более общем виде, так что они могут быть выполнены с помощью произвольной сетки. Такой подход обуславливает естественную легкость добавления новых типов сеток.

В ходе написания EMTL был реализован ряд оригинальных программных решений, подробно описываемых в рассматриваемой главе: разбиение алгоритма обновления сеток на основной цикл, поддающийся векторизации, и специализированные поправки к нему; оптимизация использования памяти с помощью упаковки используемых в разностных уравнениях коэффициентов;

увеличение параллельной эффективности путем балансировки доменов.

С целью измерения параллельной эффективности нашего кода при поддержке компании Intel были проведены тесты на системах Harpertown и Рис. 4: Слева: Ускорение на единицу моделируемого объема для различных архитектур и реализаций MPI, нормировано на время исполнения на 64 ядрах. HT Harpertown+Inniband, HT - Nehalem+Inniband (измерения в ИНТЕЛ), Xeon - кластер Т-платформы, узлы 2xXeon 5345 2.33 GHz 8 GB RAM, Inniband.

Справа: Зависимость времени исполнения задачи от числа процессоров в серии масштабируемых тестов на кластере ФМБФ МФТИ (Athlon MP 2400, Gigabit). Начальный размер задачи (сетка ) масштабируется в одном измерении пропорционально числу процессоров. Количество временных шагов равно 100.

Nehalem. В этих тестах размер сетки полагался постоянным или зависел от числа процессоров таким образом, чтобы общий объем задачи сохранялся при его увеличении. Задача представляла собой численное решение уравнений Максвелла в диэлектрике и металле. Результаты тестов приводятся на Рис. 4, где показано ускорение на единицу моделируемого объема, характеризующее масштабируемость программы. Ускорение вычисляется относительно задачи, исполняемой на 64 процессорах при размере сетки. Оба типа тестов показали линейный или суперлинейный характер масштабируемости приложения для большого числа вычислительных ядер. Также были проведены тесты на малом кластере, состоящем из 24 2-процессорных узлов с четырехядерными процессорами Xeon (2xXeon 5345 2.33 GHz 8 GB RAM).

Результаты тестов показывают, что характер производительности приложения зависит от архитектуры системы и реализации MPI.

Нами было проведено сравнение производительности с программой MEEP, разработанной в Массачусетском технологическом институте (MIT). Для этого была проведена серия масштабируемых тестов на кластере ФМБФ МФТИ (Athlon MP 2400, Gigabit). Результаты тестов показали, что производительность разработанной нами программы примерно в два раза выше, чем у MEEP (см. Рис. 4).

Третья глава посвящена изучению излучательных характеристик металлических фотонных кристаллов (ФК) и возможности создания на их основе высокоэффективных источников света. Эта глава содержит в себе литературный обзор и изложение новых результатов.

Вначале литературного обзора излагаются существующие методы расчета излучательной способности ФК. Затем рассматриваются имеющиеся в литературе результаты относительно излучательной способности металлических ФК. Отметим, что в настоящее время имеется лишь ограниченное число работ (см., напр., [9]), в которых исследуется излучательная способность ФК в видимом диапазоне. Вместе с тем, значительная часть результатов, касающихся влияния геометрических параметров ФК на спектр его излучения, получена для инфракрасного (ИК) диапазона в рамках исследований возможности применения ФК для управления спектром теплового излучения (см., напр. [3] и ссылки в этой работе). В этих работах рассматривались ФК разных геометрий: одномерные ФК, двумерные массивы рассеивателей и микрополостей, а также трехмерные ФК типа поленница. Было показано, как влияют на интенсивность излучения ФК величина периода, количество слоев ФК, оптические свойства используемого материала.

Еще раз отметим, что число работ для излучения ФК в видимом диапазоне довольно ограничено. Также не было проведено какого-либо исследования возможности использования ФК в качестве источников видимого излучения.

Этому исследованию посвящена основная часть третьей главы.

В качестве материала ФК мы используем вольфрам. Высокая селективность излучения вольфрама в видимом диапазоне выгодно отличает его от других металлов и вместе с высокой температурой плавления ( ) делает его лучшим материалом для применения в качестве источника света.

Для расчетов нами использовался метод FDTD. Для расчета излучательной способности ФК мы пользуемся способом, предполагающим расчет поглощательной способности ФК и применение закона Кирхгофа [3], связывающего ее с излучательной способностью.

Для верификации написанной программы, реализующей метод FDTD, было проведено сравнение с экспериментом, специально проведенным в компании General Electric. Экспериментальный образец представлял собой вольфрамовый ФК монослой с периодом треугольной решетки мкм. Поперечный размер рассеивающих элементов равнялся 0.2 мкм. Эксперимент проводился при комнатной температуре.

Сравнение спектра поглощения экспериментального образца с расчетами Рис. 5: Слева: Поглощение вольфрамового ФК монослоя при нормальном падении при комнатной температуре ( ). Параметры ФК: мкм, поперечный размер рассеивающих элементов 0.1 мкм. 1 - экспериментальные данные; 2 - расчет идеального ФК с учетом реальной формы элементов; 3 - расчет с учетом дефектов решетки ФК и дисперсии элементов по размерам.

FDTD показывает, что амплитуда пика поглощения в видимом диапазоне определяется степенью неупорядоченности решетки, а ИК-часть спектра поглощения определяется только размером элементов (Рис. 5). Последний факт связан с тем, что в ИК-диапазоне длина волны больше константы решетки, и поэтому электромагнитное излучение взаимодействует с ФК, как с эффективной средой.

Перейдем к результатам численных расчетов для температуры, являющейся характерной рабочей температурой вольфрамовой нити в лампах накаливания. Мы рассматриваем ФК типа гранецентрированного кубического (ГЦК) прямого опала. Такой выбор обусловлен двумя причинами.

Во-первых, прямой опал в отличие от альтернативных геометрий ФК (обратные опалы или ФК типа поленница ) обладает наилучшей термической устойчивостью при высоких температурах. Во-вторых, прямые опалы позволяют достигать малых значений фактора заполнения металлом. При этом, как нами было показано, именно при малых значениях фактора заполнения металлом достигается больший КПД светового источника.

Рассмотрим идеальный случай прямого опала, когда вольфрамовые шарики помещены в вакуум. На Рис. 5 представлена зависимость поглощательной способности двухслойного ФК от длины волны падающего света для различных значений периода решетки при сохранении величины отношения постоянной. Мы видим, что в спектре поглощения ФК присутствует ярко выраженный пик. Его положение связано с наименьшим вектором обратной решетки, величина которого равна. Это соответствует длине волны, которая определяет левый край пика.

Как видно из Рис. 5, полученные результаты в видимом диапазоне хорошо масштабируются при различных значениях и при постоянном отношении, что, на наш взгляд, связано со слабой зависимостью мнимой части диэлектрической проницаемости вольфрама в видимом диапазоне от длины волны. Таким образом, изменяя величину, можно добиться того, чтобы пик в спектре поглощения попадал в видимый диапазон.

Как показали дополнительные расчеты, при увеличении угла падения рассматриваемый пик перемещается в сторону ИК области. Увеличение радиуса шариков и числа слоев в ФК приводят к увеличению амплитуды пика, что связано с увеличением доли вольфрама в ФК. С одной стороны это приводит росту излучаемой энергии в видимом диапазоне, однако при этом падает КПД.

Нами были проведены расчеты для различных значений периода мкм. Нашей целью было нахождение таких параметров, при которых достигается максимальный КПД при сохранении величины светового потока, излучаемого через единичную площадку, не меньшим аналогичной величины для полубесконечного вольфрамового слоя. Наши расчеты показали, что искомыми параметрами ФК являются мкм, примерно равен 2.2%). При этих значениях достигается КПД,а лм/м. Отметим, что, КПД полубесконечного вольфрамового слоя равен КПД, а световой поток, излучаемый через единичную площадку, лм/м. Таким образом, применение ФК может позволить увеличить КПД примерно в 3 раза.

Дополнительные расчеты показали, что наличие неупорядоченности расположения рассеивателей в ФК и их дисперсии по размерам приводит к уменьшению пика поглощения и соответственно к уменьшению доли излучаемой энергии в видимом диапазоне. При этом сам эффект подавления излучения в ИК диапазоне сохраняется.

При изготовлении ФК образца должна использоваться внешняя матрица, фиксирующая ФК элементы. Известно лишь небольшое количество материалов, которые были бы устойчивы к высоким температурам и которые можно использовать в качестве такой матрицы (в их числе:, стабилизированный ). Однако, при нагревании эти материалы начинают поглощать, а, следовательно, и излучать в ИК диапазоне, что может привести к существенному снижению КПД. Ниже мы исследуем этот вопрос на примере Будем вначале считать, что поглощение отсутствует, и рассматривать диs= коэффициента затухания в формуле Лоренца для поглощающего (см. текст).

электрическую проницаемость при комнатной температуре ( ).

Изменение диэлектрической проницаемости внешней среды от до приводит к тому, что спектральные кривые для поглощения, полученные нами для ФК в вакууме, масштабируются на величину. Это связано с уменьшением длины волны во внешней среде в раз. Кроме того, конечность матрицы, в которую помещен ФК, приводит к появлению дополнительных осцилляций, вызванных интерференцией волн, многократно отраженных от границ между матрицей и вакуумом (Рис. 6).

При высоких температурах диэлектрическую проницаемость с хорошей точностью можно аппроксимировать с помощью формулы Лоренца где, рад/сек, рад/сек. Конкретные значения используемых параметров взяты нами из результатов экспериментов, проведенных в компании General Electric.

Для того, чтобы последовательно учесть влияние поглощения, мы полагаем рад/сек, и постепенно меняем параметр от нуля до единицы (0 соответствует отсутствию поглощения, а 1 соответствует действительной диэлектрической проницаемости ). Как видно из Рис. 6, при малых значениях спектр поглощения в основном определяется вольфрамовыми шариками. При увеличении существование матрицы вносит значительный вклад в поглощение в ИК диапазоне. При погружении ФК с найденными нами ранее оптимальными параметрами в матрицу толщины Таким образом, при росте температуры увеличение поглощения у имеющихся термически устойчивых к высоким температурам материалов приводит к полному замыванию эффекта ФК у исследованных нами образцов.

Отсюда следует, что для создания высокоэффективных источников света на основе ФК необходимо нахождение материалов со слабым поглощением при высоких температурах.

Четвертая глава посвящена изучению антиотражающих свойств тестурированных покрытий. Она состоит из литературного обзора и изложения новых результатов.

В имеющихся работах показано, что текстурированные покрытия обладают антиотражающими свойствами как в длинно-, так и в коротковолновом пределах. В первом случае это связано с тем, что оптические свойства таких покрытий подобны свойствам переходного слоя с плавно меняющимся эффективным коэффициентом преломления [10]. В другом предельном случае это связано с тем, что лучу требуются много переотражений внутри такой поверхности прежде чем окончательно отразиться обратно [11]. На примере кремниевых текстурированных покрытий было показано, что достижение низких значений коэффициента отражения возможно и в случае, когда характерный размер текстуры сравним с длиной волны [5]. Был эмпирически установлен ряд общих закономерностей касательно влияния размера и геометрии текстуры на ее антиотражающие свойства. Однако не было сделано какой-либо аналитической оценки характера убывания отражения при изменении геометрии рассеивателей. Отметим также, что в имеющихся работах, как правило, либо отдельно рассматривались длинноволновый или коротковолновый пределы, либо оптический диапазон для размеров рассеивателей порядка длины волны. Таким образом представляет интерес единое исследование текстурированных поверхностей во всем диапазоне длин волн.

В основной части главы мы приводим результаты исследования антиотражающих свойств текстурированных поверхностей во всем диапазоне размеров текстуры. Вначале мы оцениваем характер изменения отражения при изменении размеров текстуры в длинно- и коротковолновом пределах с помощью приближений эффективной среды и геометрической оптики соответственно. Затем мы исследуем антиотражающие свойства текстурированной поверхности во всем диапазоне размеров текстуры с помощью метода FDTD.

В конкретных численных расчетах мы рассматриваем подложку с нанесенными на нее стеклянными пирамидками высотой. Основаниями пирамидок являются треугольники, квадраты и круги (в этом случае пирамидка вырождается в конус). Пирамидки плотно упакованы в квадратную или треугольную решетки с периодом. При этом мы отдельно выделяем два случая:

случай полного замощения подложки основаниями пирамидок (это соответствует плотной квадратной или треугольной упаковки пирамидок с квадратными или треугольными основаниями соответственно) и случай неполного замощения (это соответствует конусам, поскольку между их круглыми основаниями всегда имеется зазор).

В длинноволновом пределе, оптические свойства текстурированной поверхности подобны свойствам промежуточного слоя с непрерывно меняющейся диэлектрической проницаемостью. Увеличение высоты пирамидок приводит к увеличению плавности изменения эффективной диэлектрической проницаемости и соответственно к уменьшению отражения.

Нами аналитически показано, что в случае полного замощения, для которого фактор заполнения (доля площади сечения пирамидок в элементарной ячейке на данной высоте) у вершин пирамидок, а у оснований коэффициент отражения при увеличении уменьшается экспоненциально (см. Рис. 7). Выводы находятся в согласии с численными расчетами, проведенными с помощью приближения эффективной среды и метода FDTD, для пирамидок с различными профилями (Рис. 7).

В случае неполного замощения подложки основаниями пирамидок между этими основаниями имеется зазор ( ). Вследствие этого в точке у функции эффективной диэлектрической проницаемости имеется разрыв. Отсюда следует, что при увеличении отражение выходит на постоянное значение, равное коэффициенту отражения между средами со значениями диэлектрической проницаемости внешней среды и текстуры (см. Рис. 7).

Перейдем теперь к рассмотрению оптических свойств текстурированной поверхности в коротковолновом пределе. На Рис. 7 представлены результаты расчетов для плотно упакованных треугольных и квадратных пирамидок (полное замощение) и конусов (неполное замощение). Результаты, полученные с помощью FDTD находятся в хорошем согласии с результатами, полученными методом трассировки лучей.

В случае полного замощения коэффициент отражения экспоненциально Рис. 7: Слева: Коэффициент отражения переходного слоя с оптическими свойствами, соответствующими (a) дифракционной решетке с интегральным профилем (TE поляризация); (b) плотно упакованным в квадратную решетку квадратным пирамидкам с (c) линейным профилем и (d) профилем (см. текст);

(e) плотно упакованным в треугольную решетку конусам. Выражения для эффективной диэлектрической проницаемости квадратных пирамидок и конусов взяты из работ [12,13].

Результаты FDTD (точки) соответствуют значениям параметров,,, Справа: Коэффициент отражения плотно упакованных пирамидок в приближении геометрической оптики. Результаты FDTD (точки) получены для.

убывает при увеличении. Мы объясняем это тем, что число переотражений внутри текстуры, необходимое для того, чтобы падающий луч поменял свое направление и отразился, линейно растет с увеличением высоты. В случае неполного замощения (например, плотно упакованные в треугольную решетку конусы на Рис. 7) при увеличении отражение вначале падает, а после достижения некоторого минимума поднимается до постоянного значения. Мы связываем это с тем, что при часть лучей после первоначального отражения от боковой поверхности движется почти параллельно ей и неминуемо попадает в зазор между основаниями, от которого эти лучи окончательно отражаются обратно.

Перейдем теперь к результатам FDTD для всего диапазона размеров текстуры. Рассмотрим вначале случай полного замощения подложки основаниями пирамидок (Рис. 8). Зафиксируем какое-либо значение и исследуем зависимость отражения от (см. соответствующие черные кривые на Рис. 8). При увеличении от до отражение падает и достигает локального минимума при, что качественно объясняется тем, эффективная диэлектрическая проницаемость в нулевом приближении не зависит от, а отражение соответствующего промежуточного слоя падает при увеличении высоты. При отражение продолжает падать, проходя ряд локальных минимумов, соответствующих тем значениям, при которых появляются очередные дифракционные порядки. Отметим, что при Reflectance, dB Рис. 8: Слева: Коэффициент отражения от плотно упакованных квадратных пирамидок при различных значениях и.

Справа: Коэффициент отражения от плотно упакованных в треугольную решетку конусов при различных значениях и.

Результаты получены с помощью метода FDTD.

, приближение эффективной среды становится неприменимым, и падение отражения обусловлено эффектами дифракции и интерференции волн, рассеянных на текстурированной поверхности. При больших значениях размах осцилляций уменьшается, и кривая выходит на предел геометрической оптики.

В случае неполного замощения, как, например, для случая плотно упакованных в треугольную решетку конусов (Рис. 8), наименьшее значение отражения достигается при.

Итак, нами установлено, что оптимальные размеры текстуры зависят от характера замощения подложки основаниями пирамидок. В случае полного замощения, когда между основаниями нет зазоров, при фиксированной величине отношения наименьшее отражение достигается при.

Это согласуется с тем, что, как мы показали, скорость убывания отражения при росте и фиксированных и в приближении геометрической оптики экспоненциальна, в то время как в приближении эффективной среды она в общем случае имеет степенной характер. При этом малая величина отражения достижима также и при порядка длины волны. В случае неполного замощения оптимальное значение имеет порядок длины волны. Наличие минимума в этом случае также находится в согласии с полученными нами оценками в приближениях эффективной среды и геометрической оптики.

В оставшейся части четвертой главы мы демонстрируем справедливость полученных выводов для кремниевых текстурированных покрытий, используемых при увеличении эффективности солнечных элементов. Мы показываем, что для достижения наименьшего отражения в видимом диапазоне наиболее эффективным решением являются нанотекстурированные покрытия.

Также нами учитывается возможная нерегулярность расположения рассеивателей, рассмотрение которой отсутствует в существующих работах. В заключение четвертой главы мы приводим реализованный нами практический метод расчета углового распределения солнечного света в течение года. С помощью этого метода рассчитан интегральный коэффициент отражения для выбранной географической точки.

3 Основные результаты и выводы работы Найдены геометрические параметры вольфрамового фотонного кристалла типа прямого опала (период решетки порядка нескольких сотен нм, фактор заполнения по вольфраму порядка нескольких процентов), погруженного в керамическую матрицу, при которых в случае отсутствия поглощения в матрице КПД источника света на основе фотонного кристалла оказывается выше в несколько раз по сравнению с обычными лампами накаливания. Показано, что у рассмотренных составляющих матрицу термоустойчивых материалов (оксид гафния) сильное поглощение при высоких температурах приводит к существенному снижению эффективности светового источника на основе фотонного кристалла, вплоть до значения КПД лампы накаливания. Сделан вывод о том, что для создания высокоэффективного источника света на основе фотонного кристалла необходим поиск материалов со слабым поглощением при высоких температурах.

Показано, что наличие дефектов в фотонных кристаллах (неупорядоченности по расположению и разброса по размерам составляющих фотонный кристалл элементов) приводит к уменьшению величины излучения в видимом диапазоне при сохранении подавления излучения в инфракрасном диапазоне. Получено количественное согласие результатов расчетов с результатами специально поставленного эксперимента. Это дает обоснование для применения использованного в расчетах метода FDTD для моделирования источников света на основе фотонных кристаллов.

С использованием трех различных методов найдены асимптотические зависимости величины отражения от геометрических параметров составляющих текстуру рассеивателей для длинно- и коротковолнового пределов. Установлено, что ключевым фактором, влияющим на оптимальный размер рассеивателей, является характер замощения подложки их основаниями. В случае полного замощения наименьшее значение отражения достигается при макроскопических размерах рассеивателей, а в случае неполного - оптимальный размер имеет порядок длины волны.

Показано, что для достижения наименьшей величины отражения в видимом диапазоне наиболее эффективным решением являются нанотекстурированные покрытия. Установлено, что наличие неупорядоченности, а также небольших зазоров в расположении составляющих эти покрытия рассеивателей не сказывается на величине отражения. При этом в случае кремниевых нанотекстурированных покрытий, используемых для увеличения эффективности солнечных элементов, практически достижимой является величина отражения порядка одного процента во всем видимом диапазоне (для нетекстурированной поверхности коэффициент отражения равен 30-40%).

Предложен новый подход программной реализации метода решения уравнений Максвелла в конечных разностях (FDTD), основанный на интегральном представлении уравнений Максвелла. В рамках этого подхода написана параллельная программа на языке программирования С++.

В ходе ее написания реализован ряд оригинальных решений: оптимизация использования памяти с помощью упаковки используемых в разностных уравнениях коэффициентов; увеличение производительности посредством разбиения алгоритма обновления сеток на основной цикл, поддающийся векторизации, и специализированные поправки к нему;

увеличение параллельной эффективности путем балансировки доменов.

Специально проведенные тесты показали, что быстродействие программы примерно в два раза выше, чем у программы MEEP, разработанной в Массачусетском технологическом институте. Также был продемонстрирован линейный характер масштабируемости программы вплоть до тысячи вычислительных ядер. Развитый при написании программы подход обеспечивает возможность расчета с одновременным использованием нескольких сеток, каждая из которых максимально приспособлена к соответствующей ей области моделируемой структуры. Это является особо важным при моделировании оптических свойств структур сложной формы с неоднородностями, много меньшими длины волны. Программа находится в открытом доступе в интернете.

В рамках метода FDTD предложен и реализован алгоритм подсеточного сглаживания для дисперсных сред, позволяющий существенно уменьшить величину численной ошибки, связанной с неточным отображением тел, обладающих произвольной формой поверхности, на прямоугольную вычислительную сетку. Его применение приводит к увеличению порядка точности метода FDTD по сравнению с обычным лестничным приближением. Это является значимым для самых разнообразных оптических приложений: при моделировании оптических свойств металлических ФК, при расчете антиотражающих свойств текстурированных покрытий и т. д.

Предложен и реализован итерационный метод для моделирования наклонного падения плоской волны на периодическую структуру в численном расчете FDTD. С его помощью можно получать результаты во всем диапазоне углов от до, а также наблюдать реальное поведение поля во времени. Это позволяет эффективно исследовать излучательные характеристики металлических фотонных кристаллов, а также антиотражающие свойства текстурированных покрытий при любом угле падения.

Основные результаты диссертации опубликованы в следующих работах:

1. I. Valuev, A. Deinega, A. Knizhnik, B. Potapkin, Creating Numerically Ecient FDTD Simulations Using Generic C++ Programming, Lecture Notes in Computer Science, 4707, 213–226 (2007).

2. A. Deinega and I. Valuev, Subpixel smoothing for conductive and dispersive media in the nite-dierence time-domain method, Optics Letters 32, 3429– 3431 (2007).

3. I. Valuev, A. Deinega, and S. Belousov, Iterative technique for analysis of periodic structures at oblique incidence in the nite-dierence time-domain method, Optics Letters 33, 1491–1493 (2008).

4. A. Deinega, I. Valuev, B. Potapkin, and Y. Lozovik, Antireective properties of pyramidally textured surfaces, Optics Letters 35, 106–108 (2010).

5. А. В. Дейнега, И. В. Конистяпина, М. В. Богданова, И. А. Валуев, Ю. Е. Лозовик, Б. В. Потапкин, Оптимизация антиотражающего слоя в солнечных батареях на основе первопринципных расчетов, Известия высших учебных заведений. Физика., 11, 13–19 (2009).

6. C. А. Белоусов, М. В. Богданова, И. А. Валуев, А. В. Дейнега, С. Л. Эйдерман, А. А. Книжник, И. Я. Полищук, Ю. Е. Лозовик, Б. В. Потапкин, Ю. А. Успенский, Э. Т. Кулатов, А. А. Титов, S. Zalyubovsky, B. Ramamurthi, Предсказательное моделирование оптических свойств металло-диэлектрических метаматериалов, Известия высших учебных заведений. Физика., 11, 20–27 (2009).

Цитируемая литература [1] Андриевский Р. А., Рагуля А. В. Наноструктурные материалы. Москва: Издательский Центр Академия, 2005.

[2] Joannopoulos J. D., G.Johnson S., Winn J. N., Meade R. D. Photonic Crystals: Molding the Flow of Light. Princeton Univ. Press, 2008.

[3] Fleming J. G., Lin S. Y., El-Kady I., Biswas R., Ho K. M. All-metallic three-dimensional photonic crystals with a large infrared bandgap // Nature. 2002. Vol. 417. Pp. 52 -55.

[4] Improved broadband and quasi-omnidirectional anti-reection properties with biomimetic silicon nanostructures / Y.-F. Huang, S. Chattopadhyay, Y.-J. Jen, C.-Y. Peng, T.-A. Liu, Y.-K. Hsu, C.-L. Pan et al. // Nature Nanotechnology. 2007. Vol. 2. Pp. 770–774.

[5] Llopis F., Tobias I. Texture prole and aspect ratio inuence on the front reectance of solar cells // J. Appl. Phys. 2006. Vol. 100. P. 124504.

[6] Taove A., Hagness S. H. Computational Electrodynamics: The Finite Dierence Time-Domain Method. Boston: Artech House, 2005.

[7] Farjadpour A., Roundy D., Rodriguez A., Ibanescu M., Bermel P. Improving accuracy by subpixel smoothing in fdtd // Optics Letters. 2006. Vol. 31. Pp. 2972–2974.

[8] Stefanou N., Yannopapas N., Modinos A. Heterostructures of photonic crystals: frequency bands and transmission coecients // Comput. Phys. Commun. 1998. Vol. 113. Pp. 49–77.

[9] Narayanaswamy A., Chen G. Thermal emission control with one-dimensional metallodielectric photonic crystal // Phys. Rev. B. 2004. Vol. 70. P. 125101.

[10] Raguin D. H., Morris G. M. Analysis of antireection-structured surfaces with continuous onedimensional surface proles // Appl. Opt. 1993. Vol. 32. Pp. 2582–2598.

[11] Sopori B. L., Pryor R. A. Design of antireection coatings for textured silicon solar cells // Sol.

Cells. 1983. Vol. 8. Pp. 249–261.

[12] Brauer R., Bryngdahl O. Design of antireection gratings with approximate and rigorous methods // Appl. Opt. 1994. Vol. 33. Pp. 7875–7882.

[13] Wu F., Whites K. W. Computation of static eective permittivity for a multiphase lattice of cylinders // Electromagnetics. 2001. Vol. 21. Pp. 97–114.

Подписан к печати 13.04. Отпечатано в отделе оперативной печати физического факультета МГУ