Содержание к диссертации
Введение
1. Методы моделирования 26
1.1. Метод прямого моделирования Монте-Карло 26
1.1.1. Общая характеристика метода 26
1.1.2. Метод и алгоритм прямого моделирования Монте-Карло в молекулярной газовой динамике 27
1.1.3. Основные процедуры и модели метода 31
1.1.4. Выбор численных параметров 39
1.1.5. Современные тенденции развития метода ПММК 41
1.2. Численный метод расчета системы уравнений Эйлера для газовой
смеси 42
2. Истечение пара от сферического источника 46
2.1. Истечение в вакуум 46
2.1.1. Постановка задачи 46
2.1.2. Результаты моделирования и их анализ 48
2.2. Истечение в фоновый газ 61
2.2.1. Постановка задачи 61
2.2.2. Результаты моделирования и их анализ 62
Заключение к главе 2 68
3. Моделирование нестационарного расширения конденсирующегося пара от сферического источника в вакуум и фоновый газ 70
3.1. Нестационарное расширение конденсирующегося пара от сферического источника в вакуум (модель ПММК) 70
3.1.1. Математическая модель 70
3.1.2. Результаты моделирования и их анализ 75
3.2. Нестационарное расширение конденсирующегося пара от сферического источника в вакуум (приблтісенная газодинамическая модель) 83
3.2.1. Математическая модель 83
3.2.1.1. Основные допущения 83
3.2.1.2. Модели констант скоростей кинетических процессов 84
3.2.1.3. Модель энергообмена 88
3.2.1.4. Основные уравнения 89 3.2.2. Результаты моделирования и их анализ 90
3.3. Нестационарное расширение конденсирующегося пара от
сферического источника в фоновый газ (модель ПММК) 98
3.3.1. Математическая модель 98
3.3.2. Результаты моделирования и их анализ 99
3.4. Моделирование процессов образования и роста кластеров в
быстро охлажденной однородной смеси пара и инертного газа 106
3.4.1. Схема, параметры и основные стадии процесса получения нанокластеров кремния методом ИЛА частиц аэрозоля 106
3.4.2. Кинетическая модель процесса образования нанокластеров кремния (пространственно-однородная релаксация) 109
3.4.3. Результаты моделирования процессов образования нанокластеров Si в рабо чей камере 114
Заключение к главе 3 117
4. Моделирование газодинамических процессов ИЛА при переходе к режиму взрывного вскипания 120
4.1. Математическая модель 121
4.1.1. Тепловая модель мишени 121
4.1.2. Тепловая модель оторвавшегося слоя 123
4.1.3. Модель движения пара 124
4.2. Результаты моделирования и их анализ 125
4.2.1. Оптическая модель слоя 1 125
4.2.2. Оптическая модель слоя 2 129
4.3. Обсуждение результатов 130
Заключение к главе 4 131
Заключение 132
Литература
- Метод и алгоритм прямого моделирования Монте-Карло в молекулярной газовой динамике
- Результаты моделирования и их анализ
- Нестационарное расширение конденсирующегося пара от сферического источника в вакуум (приблтісенная газодинамическая модель)
- Тепловая модель оторвавшегося слоя
Введение к работе
Диссертация посвящена математическому моделированию газодинамических и физических процессов, сопровождающих абляцию твердого вещества (на примере кремния) в вакууме и в фоновом газе импульсами лазерного излучения (ЛИ) наносекундного диапазона умеренной интенсивности. Явление импульсной лазерной абляции (ИЛА) представляет собой стимулированный лазерным воздействием унос вещества с облучаемой поверхности, сопровождающийся образованием парогазового облака и его последующим разлетом в окружающее пространство. Выбор кремния в качестве исследуемого вещества обусловлен тем, что он является основным элементом электроники и представляет особый интерес для нанотехнологий.
Выбор темы диссертации стимулирован широким использованием импульсной лазерной абляции в ряде современных технологий, в том числе технологий получения нанокластеров и наноматериалов. Решающее значение для выбора конкретных задач данной работы имели исследования процессов получения нанокластеров кремния при наносекундной ИЛА мелких аэрозольных частиц, выполненные Центром перспективных исследований СПбГПУ по Соглашению с Samsung Electronics Co., Ltd. в 2004 году [1] и продолженные в последующие годы в рамках гранта INTAS №03-51-5208 "Механизмы образования кластеров при импульсной лазерной абляции" (2004-2007) и плана фундаментальных исследований ЦПИ СПбГПУ.
Импульсная лазерная абляция и ее приложения. Первыми исторически возникли такие области применения ИЛА как сварка, резка, сверление, обжиг, упрочнение твердых поверхностей [2-5]. В 70-е годы в связи с идеей лазерного термоядерного синтеза появился интерес к облучению материалов мощными наносекундными импульсами, что повлекло за собой большое количество экспериментальных работ [например, 6-8]. Новый всплеск интереса к наносекундной лазерной абляции возник в середине 80-х годов в связи с открытием высокотемпературных сверхпроводников [9], развитием технологий напыления тонких пленок [10, 11], разработкой эффективного источника кластерных пучков на основе лазерной абляции [12, 13] и обнаружением фуллеренов в лазерной плазме [14]. В настоящее время наносекундная лазерная абляция является наиболее пригодным методом напыления тонких пленок со сложной стехиометрией и микрообработки хрупких материалов [10, 11, 15] и мощным инструментом изготовления компонентов микроприборов и микрооптических устройств, очистки поверхностей, сверления отверстий микронного размера и др. [15]. Наносекундные лазерные импульсы все шире используются в медицине (офтальмология, дерматология, клеточная хирургия) [16].
Метод и алгоритм прямого моделирования Монте-Карло в молекулярной газовой динамике
Метод Монте-Карло обычно определяют как метод моделирования случайных величин с целью вычисления характеристик их распределения. Метод Монте-Карло, получивший в настоящее время применение в самых различных областях науки, был развит группой физиков и математиков, работавших в Лос-Аламосе, включая Джона фон Неймана и Стенли Улама. Датой рождения метода принято считать 1949 год, когда появилась статья под названием "The Monte Carlo method", в которой метод Монте-Карло впервые излагался систематическим образом [155].
Первоначально метод Монте-Карло использовался в основном для решения задач нейтронной физики, где традиционные численные методы оказались малопригодными. Далее его влияние распространилось на широкий класс задач статистической физики, теории массового обслуживания, теории игр, математической экономики, передачи сообщений при наличии помех, экологии и ряд других. Весьма характерным является прогресс метода Монте-Карло в приложении к задачам динамики разреженного газа. В настоящее время метод Монте-Карло практически вытеснил все иные подходы к решению задач в этой области.
В настоящее время теория методов Монте-Карло является вполне сформировавшейся научной дисциплиной. Имеется ряд прекрасных монографий и учебных пособий, в которых рассмотрены общие вопросы метода, в том числе его математическое обоснование, а также примеры разнообразных приложений метода, там же можно найти обширную библиографию [156-160].
Метод прямого моделирования Монте-Карло является основным методом молекулярной газовой динамики и в настоящее время широко используется для решения самых разнообразных задач динамики разреженного газа, включая задачи аэрокосмической техники, динамики верхней атмосферы планет, вакуумной техники и т.д. Прекрасное изложение метода содержится в монографиях Г. Бер да [93, 161].
Метод ПММК является техникой моделирования на ЭВМ реального течения газа достаточно большим числом моделирующих молекул. В памяти компьютера запоминаются координаты и скорости моделирующих молекул, изменение которых в течение времени обусловлено межмолекулярными столкновениями и взаимодействием с границами физического пространства.
Моделируемый объем физического пространства разбивается на ячейки. Такими ячейками могут быть либо малые объемы с заданными границами, либо просто множество точек. В последнем случае, молекула находится в данной ячейке, если она ближе всего к точке, определяющей эту ячейку. Размеры ячеек должны быть такими, чтобы изменение параметров течения в каждой ячейке было малым.
Метод ПММК использует принцип расщепления непрерывного процесса движения и столкновений молекул в разреженном газе на два последовательных независимых этапа на каждом малом временном шаге At, величина которого во всей области должна удовлетворять условию где xav - среднее время между столкновениями частиц.
Моделирование начинается с установления начального состояния моделирующих молекул внутри области и задания граничных условий. Моделирование протекает во времени (по временным шагам) и всегда воспроизводит некоторый нестационарный процесс. Стационарное решение, если оно существует, достигается при заданном начальном состоянии и стационарных граничных условиях в результате процесса установления за некоторое время tst. При решении стационарной задачи методом установления выбор начального состояния, вообще говоря, является произвольным.
Конечная величина памяти и быстродействия ЭВМ определяют ограничения на общее число моделирующих молекул и количество ячеек. Макроскопические параметры газа, определяемые параметрами моделирующих частиц в ячейках на данном временном шаге, представляют собой результат выполненного статистического испытания (процедуры моделирования). Ввиду относительно небольшого числа моделирующих частиц в ячейке статистические флуктуации средних параметров в ячейках на временных шагах весьма велики. Поэтому при решении стационарных задач после установления течения для уменьшения статистической погрешности определения параметров потока увеличивают время осреднения результатов моделирования (увеличивают объем выборки), выбирая при t tst для осреднения достаточно большой интервал времени Atav » At. Этот интервал может состоять из необходимого числа последовательных временных шагов, либо из временных шагов, выбираемых с некоторым интервалом (например, с интервалом через 5-Ю шагов). В последнем случае достигается статистическая независимость результатов, получаемых на используемых для осреднения временных шагах [93]. Статистическая ошибка в результатах обратно пропорциональна корню квадратному из числа испытаний, поэтому достижение необходимой точности моделирования требует осреднения по достаточно большому числу временных шагов Q. Величина Q зависит от числа моделирующих частиц, параметров задачи и целей исследования.
Результаты моделирования и их анализ
Одним из основных направлений развития метода прямого статистического моделирования является его распространение на исследование течений разреженного газа, сопровождающихся разнообразными химическими реакциями. В том числе в рамках метода ПММК активно разрабатываются модели процессов образования кластеров [140-145] (обзор данных работ представлен во Введении). Отметим также, что одной из современных тенденций в области моделирования течений при импульсной лазерной абляции твердых материалов с учетом процессов конденсации в потоке является развитие комбинированных методов молекулярной динамики и ПММК. Так, в работе [143] предложена комбинированная модель, в которой использованы методы молекулярной динамики при моделировании взаимодействия лазерного излучения с веществом и выброса атомов и кластеров с поверхности и метод ПММК для описания расширения эжектированного облака с учетом процессов слипания и развала кластеров.
Важным элементом моделей кластерообразования в рамках ПММК является описание энергообмена между поступательными и внутренними степенями свободы многоатомных частиц. В большинстве исследований используется модель Ларсена-Боргнакке [165]. В последнее время появляются различные модификации данной схемы, такие, например, как модель с дискретными внутренними энергиями [166] и "vibrationally favored" дискретная модель [167], которые дают более точные значения скоростей химических реакций при низких температурах.
Ввиду высокой трудоемкости метода ПММК при расчете околоконтинуальных течений в последнее время активно разрабатывается новый класс более экономичных схем "time relaxed Monte Carlo" (TRMC) [168, 169]. TRMC схемы основаны на формальном распространении решения уравнения Больцмана, содержащего только положительные члены, и его последующей вероятностной интерпретации. В работе [168] предложена и протестирована TRMC схема для случая пространственно-однородной релаксации, а в [169] - для плоского течения Куэтта.
Одним из достоинств метода ПММК является простота его распространения на параллельные вычисления. К настоящему времени имеется уже достаточно большое количество работ, посвященных использованию параллельных алгоритмов при решении задач динамики разреженного газа [170, 171].
Для расширения возможностей использования метода прямого моделирования Монте-Карло в различных инженерных приложениях Бердом и др. ведутся разработки по созданию пользовательского программного пакета, предназначенного для решения широкого класса задач динамики разреженного газа, в том числе многомерных задач со сложной геометрией [172].
Численный метод расчета системы уравнений Эйлера для газовой смеси
В разделе 3.2 диссертационной работы рассматривается задача о нестационарном сферическом разлете пара в вакуум от звукового источника при наличии процессов кластерообразования на основе математической модели, построенной в рамках газодинамического приближения Эйлера.
Математическая модель течения включает уравнения баланса массы для пара в целом и отдельных компонент (атомов и кластеров), изменения количества движения, энергии и состояния. Учитывая, что энтальпия смеси h = Y hkkck = cpT + hkkck (hk=cpkT + hk - энтальпия кластера к, ск к к объемная доля кластера к), систему уравнений можно привести к следующему безразмерному виду (см. раздел 3.2): Здесь r - радиальная координата, p - плотность, p - давление, и - скорость, у = с I cv - отношение удельных теплоємкостей для смеси атомов и кластеров, Wk - безразмерная скорость образования массы к-то компонента в единице объема {к = \,...,кт; кт - максимальный размер кластера в смеси). Характерной особенностью рассматриваемого течения при наличии процессов конденсации в потоке является появление у поверхности источника узкой перегретой дозвуковой зоны, что делает невозможным использование маршевых методов, применяемых в случае разлета простого нереагирующего газа. dF_ dU
Ввиду этой особенности для численных расчетов была использована явная конечно-разностная схема расщепления векторов потоков первого порядка точности [173]. Данный метод предполагает, что поток F обладает и свойством однородности, то есть одновременно выполняется А =
F - AU. Матрицу А можно привести к диагональному виду: A=RDL, где R -матрица правых собственных векторов, L - матрица левых собственных векторов, D - диагональная матрица, элементами которой являются собственные числа матрицы А. Пусть матрица D+ получена из матрицы D заменой отрицательных элементов нулями, а матрица D - заменой положительных элементов нулями. Тогда поток F можно представить в виде: F = F+ + F , F±=RD±LU. (1.2.3) Для системы уравнений (1.2.1) получаем следующее выражение для векторов потоков F1: Аппроксимация уравнений (1.2.1) выглядит следующим образом (1.2.5) Выбор численных параметров
Выбираемыми параметрами данной численной схемы являются безразмерные шаг по времени Ai = Atl{Rwlи,) и шаг по пространству Ar = Ar I Rw (Rw - радиус источника, и. - скорость звука в атомарном паре на поверхности источника). Условие устойчивости явной схемы накладывает \Х \At следующее ограничение на временной шаг: CFL = LJ! — 1 (в нашем Аг случае максимальное собственное число матрицы А равно (и+с)). Для задачи об истечении простого неконденсирующегося пара это ограничение является определяющим. Однако мы рассматриваем случай, когда в потоке протекают реакции кластерообразования, а именно: трехчастичная атомная рекомбинация, диссоциация димеров, слипание кластеров с атомами при парных столкновениях и испарение кластерами атомов. Причем характерное время кинетических процессов кластерообразования существенно меньше (К конвективного времени. Так, число Данкеллера Da = —— - = -= р. К хш характерная скорость образования компонента) для реакции трехчастичной атомной рекомбинации для варианта RW=10 6M составляет 1.6 103 (во всех расчетах температура на поверхности источника полагалась равной Г =6000К). В связи с этим на шаг по времени накладывается еще одно ограничение, которое и является определяющим в рассматриваемых условиях: Дт тКин, где ткин - характерное время изменения состава газовой смеси за счет процессов роста и развала кластеров. Таким образом был определен максимально возможный временной шаг для каждого из вариантов расчета.
Нестационарное расширение конденсирующегося пара от сферического источника в вакуум (приблтісенная газодинамическая модель)
В данном разделе описывается модель нестационарного истечения пара в вакуум от сферического источника при наличии процессов образования кластеров, основанная на приближении Эйлера.
Целесообразность разработки наряду с моделью ПММК приближенной газодинамической модели определяется следующими соображениями:
1. Газодинамическая модель значительно более эффективна (менее трудоемка) в вычислительном отношении.
2. Газодинамическая модель, в отличие от модели ПММК, дает явное описание скоростей основных физических процессов и процессов энергообмена. Это важно при использовании модели для оценки относительной роли отдельных процессов и механизмов энергообмена.
На пути создания качественной газодинамической модели для рассматриваемого типа течения имеются, однако, существенные трудности. Рассматриваемый диапазон условий требует, вообще говоря, использования модели вязкого теплопроводного газа для смеси, состоящей из атомов и кластеров. В интервале температур ТКЗОО-бОООК данные о коэффициентах переноса для пара кремния отсутствуют даже для равновесных условий. Почти столь же неопределенна ситуация в отношении моделей констант скоростей процессов атомной рекомбинации, диссоциации, слипания частиц и испарения кластеров, а также процессов энергообмена между различными степенями свободы кластеров. Все эти обстоятельства определяют существенно приближенный характер описываемой ниже газодинамической модели. В части описания процессов образования кластеров модель является "предельной" и соответствует случаю разлета кластеров с постоянной колебательной энергией ("горячих" кластеров).
Математическая модель 3.2.1.1. Основные допущения 1. Рассматривается односкоростная модель течения. Это допущение обосновывается результатами, приведенными в п. 3.1. 2. Пар состоит из атомов и кластеров. 3. Пар является невязким и нетеплопроводным. 4. Для описания столкновений атомов и кластеров используется модель твердых сфер. 5. Параметры течения на поверхности сферы соответствуют числу Маха М=1, определенному по "замороженной" скорости звука. 6. Пар является разреженным газом, что позволяет использовать уравнение состояния идеального газа. 7. Рассматривается диапазон условий, когда начальный состав пара близок к атомарному (массовые доли димеров и более крупных кластеров малы). 8. Рассматривается диапазон условий, при которых массовая доля атомов в паре существенно превышает суммарную массовую долю кластеров в процессе расширения. 9. Кинетическая модель учитывает следующие процессы: трехчастичную атомную рекомбинацию и диссоциацию димеров AM&Ai -» A2+Aj, (3.2.1) где А І - атом, А2 - димер; процессы слипания при парных столкновениях кластеров с атомами Ак+А, Ам,к 2; (3.2.2) процессы испарения (диссоциации) кластеров к Ъ Ак Ак.,+ А,. (3.2.3) 10. Используется приближенная модель энергообмена, описываемая ниже в п. 3.2.1.3.
Для описания процесса рекомбинации (3.2.1) используем модель трехчастичной атомной рекомбинации на базе упрощенной модели тройных столкновений [136]. Основными допущениями данной модели являются: 1. Модель предполагает, что реакция происходит при каждом тройном столкновении (3.2.1), которое трактуется как совокупность двух последовательных двойных столкновений - образования квазимолекулы Аг и ее столкновения с третьей частицей В. Время жизни квазимолекулы составляет около 10" с, что характерно для сталкивающихся молекул. 2. Столкновения рассматриваются в рамках модели твердых сфер. 3. Массы реагирующих частиц, а также радиусы их взаимодействия различаются незначительно.
Тепловая модель оторвавшегося слоя
На рис. 3.2.4 приведены зависимости dck I dt от г в различные моменты времени г для различных процессов, входящих в систему кинетических уравнений, для варианта /?W=10 8M. ЭТОТ вариант по параметрам источника является наиболее разреженным (Кп=1) из всех рассмотренных вариантов. Представленные данные позволяют сделать некоторые выводы: Процесс образования димеров (на стадии быстрого роста с2 с последующей стабилизацией значения с2) определяется практически только процессом трехчастичной атомной рекомбинации (Кг). Для более плотных вариантов расчета этот вывод должен быть еще более уверенным. Роль процесса диссоциации (Kd) димеров в данных условиях пренебрежимо мала. Процесс образования тримеров в зоне, прилежащей к источнику, определяется процессом А2+А -» А3. С удалением от источника на некотором расстоянии от него достигается баланс между процессами А2+А - А3 (Ksl2) и А3 - А2+А (Kei). Далее процесс испарения становится определяющим, что и приводит к уменьшению съ. Аналогичная картина имеет место и для кластеров к=А. Сравнение результатов, полученных с помощью модели ПММК и газодинамической модели Оценка правомочности принятых допущений
Расчеты по модели ПММК (раздел 3.1) показали, что внутренняя энергия кластеров убывает с удалением от поверхности источника, кластеры покидают расчетную область достаточно холодными (рис. 3.1.5). В данном разделе в рамках приближения Эйлера рассмотрена предельная модель "горячих" кластеров (колебательная энергия кластеров полагается неизменной во времени и в пространстве), что приводит к переоценке роли процессов испарения кластеров и ограничивает их рост.
Сделанное в данной приближенной газодинамической модели предположение о том, что течение является односкоростным, полностью подтверждается расчетами по методу ПММК (рис. 3.1.4 а). Сравнение численных результатов
Анализ профилей газодинамических параметров, полученных по двум рассмотренным моделям расширения конденсирующегося пара в вакуум, показал, что за исключением приповерхностной зоны в несколько длин свободного пробега (кнудсеновского слоя) обе модели демонстрируют качественно одинаковый характер влияния процессов конденсации на газодинамическую картину разлета пара (температура и скорость потока при наличии процессов образования кластеров выше, чем в случае расширения атомарного пара, а плотность ниже). Таким образом, рассмотренная в данном разделе приближенная модель качественно верно описывает влияние кинетических процессов кластерообразования на газодинамику истечения от сферического источника в вакуум.
Сравним теперь полученные с помощью обеих моделей данные по распределениям объемных долей кластеров для варианта RW=\0 6M, 7 6000К при т=5 (рис. 3.2.2 d, 3.2.3 b и 3.2.5). В расчетах по приближенной газодинамической модели объемная доля димеров "замораживается", а профили с для более крупных кластеров резко убывают с удалением от поверхности источника. В случае же ПММК модели эффект "замораживания" имеет место для всех кластеров. Как отмечалось выше, такое поведение с при к 3 для приближенной газодинамической модели объясняется высокими колебательными температурами кластеров, поддерживающими высокие скорости процессов испарения. димеров в кластеры прекращается, оказывается завышенной
Так как с отходом от поверхности источника трансформация более крупные практически их доля несколько по сравнению с результатами по ПММК модели. Значения же объемных долей кластеров к 3 вблизи источника для обеих моделей приближено совпадают. к = 2 Основные выводы.
1. В рамках газодинамического приближения Эйлера разработана математическая модель нестационарного расширения пара от внезапно включенного сферического источника с учетом процессов образования кластеров. Модель учитывает процессы трехчастичной атомной рекомбинации и диссоциации, процессы слипания частиц при парных столкновениях атомов с кластерами и процессы испарения кластерами атомов. Модель предполагает приближенное описание энергоообмена в процессе образования и распада кластеров: колебательная температура кластеров считается постоянной. Модель является "предельной" и позволяет оценить ход процессов образования и распада кластеров в случае, когда кластеры при разлете остаются "горячими".
2. Проведены численные исследования для случая разлета пара кремния от источника с радиусом /?W=10"8-10"5M при Г бОООК и т 10.
