Моделирование нагрева полупрозрачного материала импульсным лазером

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

Нагрев кремниевой пластины лазерным лучом

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

Импульсный лазерный нагрев двух слоёв полупрозрачного материала на непрозрачной пластине.

Поскольку лазерный луч направлен по нормали к плоской поверхности пластины, падающее излучение будет распространяться в направлении, параллельном падающему лучу. На границе раздела материалов луч будет частично отражаться, но преломления или дифракции не будет. Толщины обоих слоёв намного превышают длину волны лазера, поэтому можно предположить, что длина когерентности много меньше толщины слоя. Поставленную задачу можно решить, если воспользоваться законом Бугера — Ламберта — Бера, который описывает ослабление луча в полупрозрачной среде. Для решения соответствующего уравнения в среде COMSOL Multiphysics^® предназначен интерфейс Radiative Beam in Absorbing Media. Однако, поскольку в системе имеется отражение, появляются некоторые особенности, которые нужно будет учесть.

Анализ физических процессов и настройка расчётной модели

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

Для начала давайте посмотрим, как лазерный луч попадает от источника на пластину, проходя через пустое пространство вдоль оси z. Допустим, что мощность источника составляет 40 Вт, длина волны лазера — 700 нм, а плотность энергии в луче описывается гауссовым распределением со стандартным отклонением 1.5 мм. Лазер работает в течение 75 мс, а затем отключается на 25 мс, то есть нагрев осуществляется в импульсном режиме с периодом 100 мс. Подобный ступенчатый характер изменения тепловой нагрузки во времени описывается с помощью интерфейса Events, в котором через узел Discrete State задаётся переменная состояния ONOFF, принимающая дискретные значения 0 или 1.

Мы не будем явно моделировать ни лазерный источник, ни движение луча через свободное пространство; проанализируем только взаимодействие излучения с веществом. На границе верхнего слоя материала (показатель преломления n_{top}=2.4) луч будет частично отражаться. Этот эффект обусловлен резким изменением показателя преломления на границе раздела двух сред и описывается уравненем Френеля:

Хотя это уравнение записано для комплексных показателей преломления, в наших расчётах целесообразно рассматривать только вещественную часть показателя преломления, поскольку мнимая часть крайне мала. Если предположить, что на границе раздела нет никакого поглощения (например, по причине малой толщины слоя поглощающего материала), тогда коэффициент пропускания будет равен T=1-R. Теперь мы располагаем всей информацией, которая необходима для настройки узла Incident Intensity интерфейса Radiative Beam in Absorbing Media, как показано на скриншоте ниже.

Окно настройки узла Incident Intensity.

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

где k — мнимая часть показателя преломления, а \lambda_0 — длина волны лазера в свободном пространстве. Коэффициент поглощения может зависеть от температуры, но мы будем считать его постоянной величиной. Зная распределение плотности энергии луча по верхней поверхности, можно рассчитать плотность энергии луча во всём домене.

На границе раздела между верхним и нижним слоями луч снова будет частично отражаться в соответствии с уравнением Френеля. Отраженное излучение можно смоделировать с помощью уже имеющегося в модели физического интерфейса Radiative Beam in Absorbing Media, просто добавив ещё один узел Incident Intensity. В этом интерфейсе можно использовать любое количество узлов Incident Intensity. С каждым таким узлом связана дополнительная встроенная переменная rbam.I1, rbam.I2, ..., и так далее. В окне настройки второго узла Incident Intensity мы можем задать пользовательское соотношение для распределения плотности энергии в луче, зная плотность энергии первого луча и коэффициент отражения Френеля. Меняя знак вектора направления луча, можно учесть частичное отражение излучения на границе раздела, как показано на скриншоте ниже. Теоретически, на верхней границе произойдёт повторное отражение этого луча, но поскольку плотность энергии повторно отражённого излучения будет мала, мы можем им пренебречь.

Окно настройки второго узла Incident Intensity , используемого для учёта отражения от границы раздела диэлектриков.

Теперь проанализируем, как луч проходит второй слой полупрозрачного материала. Поскольку плотность энергии при прохождении через границу раздела меняется, мы должны добавить второй интерфейс Radiative Beam in Absorbing Media и задать интенсивность падающего излучения, зная значение коэффициента пропускания Френеля и интенсивность первого луча из первого интерфейса Radiative Beam in Absorbing Media.

Окно настройки узла Incident Intensity интерфейса Radiative Beam in Absorbing Media , с помощью которого рассчитывается излучение в нижнем домене.

Наконец, давайте посмотрим, что происходит, когда излучение достигает нижней части второго слоя и попадает на кремниевую подложку. Будем считать, что кремниевая пластина легирована и обладает высокой поглощающей способностью. Поскольку всё излучение, достигающее этой границы, будет поглощено на достаточно малом расстоянии, можно утверждать, что излучение поглощается на границе. Для описания этого эффекта можно использовать граничное условие непрозрачной поверхности Opaque Surface. Это завершающая настройка модели переноса лазерного излучения через вещество. Использованная нами комбинация интерфейсов полностью описывает прохождение лазерного луча через материалы модели. Теперь перейдём к настройке тепловой модели.

Моделирование нестационарного нагрева

В начальный момент времени температуры пластины равна 300 К. Теплота во всех доменах переносится теплопроводностью. Мы будем считать, что контактное термическое сопротивление на границах раздела материалов отсутствует, то есть поле температуры и плотность теплового потока являются непрерывными. В COMSOL Multiphysics данное предположение используется по умолчанию, однако в случае, когда оно неприменимо, можно использовать узлы для описания тонкого теплопроводного слоя Thin Layer или контактного термического сопротивления Thermal Contact.

Толщина слоёв 100 мкм достаточно велика, поэтому распределение температуры можно описать классическим уравнением Фурье, при этом отметим, что теплопередача в наноструктурах является активной областью исследований среди пользователей COMSOL; см., например, статью приглашённого автора в блоге COMSOL “Hydrodynamic Thermal Transport in the Kinetic-Collective Model”.

Для постановки граничных условий в тепловой задаче примем, что пластина лежит на идеально теплоизолирующем основании в вакуумной камере. Это означает, что конвективное охлаждение пластины отсутствует, но при этом имеет место теплообмен излучением между пластиной и стенками камеры, которые имеют постоянную температуру 300 К. Также предположим, что температура пластины увеличится всего на несколько сотен градусов, поэтому длина волны теплового излучения будет намного больше длины волны лазера. Это предположение позволяет нам использовать двухдиапазонную модель теплообмена излучением. Падающее излучение лазера полностью описывается интерфейсами Radiative Beam in Absorbing Media. Тепловое излучение пластины в длинноволновом диапазоне (обусловленное увеличением температуры пластины относительно температуры стенок камеры) можно смоделировать с помощью интерфейса Surface-to-Surface Radiation, объединённого с интерфейсом Heat Transfer in Solids. Интерфейс Surface-to-Surface Radiation позволяет рассчитать угловые коэффициенты излучения для всех поверхностей, участвующих в теплообмене.

Стоит отметить, что в данном случае теплообмен излучением между поверхностями реализуется только вблизи небольшого внутреннего угла, где слои выступают над пластиной; на всех остальных участках расчётной области угловые коэффициенты излучения в окружающую среду равны единице. Чтобы немного упростить задачу, можно было бы отказаться от интерфейса Surface-to-Surface Radiation, заменив его граничным условием Surface-to-Ambient Radiation в интерфейсе Heat Transfer in Solids. Однако, поскольку отличия в расчётном времени и в результатах моделирования незначительны, мы воспользуемся более универсальным и точным методом расчёта угловых коэффициентов излучения с помощью интерфейса Surface-to-Surface Radiation.

Особое внимание необходимо уделить построению расчётной сетки. В интерфейсах Radiative Beam in Absorbing Media решаются дифференциальные уравнения первого порядка, и по умолчанию используется первый порядок дискретизации искомых переменных. Поскольку коэффициенты поглощения нам известны, мы знаем, что плотность энергии будет значительно меняться по толщине двух слоёв. Также известно, что изменение плотности энергии в пределах луча довольно плавное. При этих условиях мы можем в пределах слоёв использовать структурированную сетку, с прямоугольными вытянутыми элементами. Конечно, по мере нарастания сложности задачи всегда следует анализировать зависимость решения от параметров сетки и допуска решателя, о чём мы рассказывали в предыдущей статье блога “Введение в моделирование нестационарной теплопередачи в твёрдых телах с помощью COMSOL Multiphysics^®”.

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

Распределение температуры пластины вдоль вертикальной оси z.

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

Изменение температуры при постоянном и зависящем от температуры коэффициенте поглощения.

Заключительные замечания

В этой статье мы описали алгоритм решения задачи о нагреве полупрозрачного материала. Лазер, являющийся в этой задаче источником теплоты, моделируется с помощью интерфейсов Radiative Beam in Absorbing Media, которые позволяют учесть полупрозрачность материала для лазерного излучения, а также отражения на границе раздела сред. Импульсный источник тепла моделируется с помощью интерфейса Events, а длинноволновое инфракрасное излучение описывается с помощью интерфейса Surface-to-Surface Radiation. Данный подход применим к задачам из области обработки полупроводников или любым другим задачам, связанным с прохождением излучения через полупрозрачные материалы.

Если вы интересуетесь решением подобных задач, пожалуйста, скачайте пример описанной в этой статье модели по представленной ниже ссылке: