автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Математическое моделирование импульсного лазерного воздействия на металлы

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

ЛОБОК Максим Геннадьевич

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ИМПУЛЬСНОГО ЛАЗЕРНОГО ВОЗДЕЙСТВИЯ НА МЕТАЛЛЫ

Специальность 05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

7 НОЯ 2013

Москва - 2013

005537181

Работа выполнена в федеральном государственном бюджетном учреждении науки «Институт прикладной математики имени М.В. Келдыша РАН»

Научный руководитель:

Официальные оппоненты:

Ведущая организация:

доктор физико-математических профессор, заведующий сектором Мажукин Владимир Иванович

наук,

Галанин Михаил Павлович, доктор физико-математических наук, профессор, федеральное государственное бюджетное учреждение науки «Институт прикладной математики имени М.В. Келдыша РАН», заведующий отделом

Каптильный Александр Григорьевич,

кандидат физико-математических наук, федеральное государственное бюджетное учреждение науки Объединенный институт высоких температур РАН, старший научный сотрудник

Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования

«Национальный исследовательский университет «МЭИ»

Защита диссертации состоится «26» ноября 2013 г. В 11-00 часов на заседании диссертационного совета Д 212.141.15 при Московском государственном техническом университете имени Н.Э. Баумана по адресу: Москва, Рубцовская наб., д. 2/18, ауд. 1006 л.

С диссертацией можно ознакомится в библиотеке Московского государственного технического университета имени Н.Э. Баумана.

Автореферат разослан «¿3 г.

Ученый секретарь диссертационного совета, кандидат технических наук, доцент

А.В. Аттетков

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

Актуальность работы. Большинство технологических приложений импульсной лазерной обработки материалов связано с началом фазовых превращений 1-го рода: плавления-затвердевания и испарения-конденсации. Соотношение глубин расплава и испарённого слоя играет важную роль в процессе микрообработки материалов.

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

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

Задачи, подлежащие исследованию:

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

2. Разработка вычислительного алгоритма для численного решения задач Стефана, сопряженных с задачами радиационной газовой динамики (РГД), с явным выделением подвижных межфазных границ и фронтов ударных волн.

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

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

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

Научная новизна. В диссертации получены следующие новые научные результаты, выносимые на защиту:

1. Методами математического моделирования исследовано нелинейное взаимодействие фазовых фронтов плавления и испарения, оказывающего количественное и качественное влияние на процессы в зоне облучения.

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

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

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

В милли - микросекундном диапазоне длительности перераспределение

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

Апробация результатов работы. Основные результаты диссертационной работы обсуждались и докладывались на следующих 10-и конференциях: Международная научная конференция «Ломоносов - 2004» (Москва, 2004), I European Summer School (Saint - Etienne, 2004), III European Summer School (Saint - Etienne, 2006), E-MRS 2006 Spring Meeting, (Nice, 2006), IV Международный научный семинар «Математические модели и моделирование в лазерно-плазменных процессах» (Москва, 2007), Third International Conference Computational methods in applied mathematics (Minsk, 2007), III International Conference on Adaptive Modeling and Simulation ADMOS 2007 (Göteborg, 2007), E-MRS 2008 Spring Meeting (Strasbourg, 2008), VI International Conference on Photo-Excited Processes and Applications (Sapporo, 2008), International Conference «Advanced Laser Technologies» (Siofok, 2008).

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

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

Структура и объем диссертации. Диссертация состоит из введения и трех глав, заключения и списка литературы. Диссертация изложена на ИЗ страницах, включает библиографический список из 146 наименований работ, иллюстрирована 20 рисунками.

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

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

Математическое описание процессов лазерного нагрева, плавления и испарения металлов производится в рамках совмещенного варианта задачи Стефана (V.l. Mazhukin, A.A. Samarskii; 1994), объединяющего классическую и однофазную постановки

W—Н(Т) = рС„(Т)Т, *-{,./}; (1)

Ao<*<r,,(OUr(((r)<.v<rh(f), т > 0; Т(х,0)-Т0,

здесь x,t - пространственная и временная переменные, Т- температура, IV - тепловой поток, //('/'), Ср('/'), А(7')- энтальпия, теплоемкость и теплопроводность конденсированной среды, £ = {*,(} - индексы, означающие принадлежность величин к твердой и жидкой фазам соответственно. Левая граница х-х0 мишени полагалась теплоизолированной VV(r,.t„) = 0. На границе раздела -т=Гл, твёрдой и жидкой фаз в классическом варианте задачи Стефана предполагается выполнение эмпирического условия равенства температуры раздела фаз Г,. равновесной температуре плавления Г, и так называемого дифференциального условия:

.<-= Г, ,(f): W,t=W,-Wt = Lmp,vit-, T=T.=T„r (2)

Процесс поверхностного испарения описывается в приближении однофазного варианта задачи Стефана с учетом кинетики испарения на кнудсеновском слое. На границе раздела фаз Г,, выписываются три закона сохранения и два соотношения на кнудсеновском слое: •г = Г4„('): -С,„ + <:г'С;

Те-Т1.(Т1,..му,р1.-р,.(р„,М). 4

дН(Г)_ д\У а/ дх

Число Маха М в общем случае определяется из решения уравнений газовой динамики при М = 1 (D. Crout; 1936)

7^. -С, *rt„, = С, *рн, С, = 0.633,С, =0.326.

Математической особенностью рассматриваемой модели является наличие двух подвижных фазовых границ Г,,(f) и ГЛ(/), которые соответствуют поверхностям раздела фаз. Положение фазовых границ заранее неизвестно и определяется в ходе численного решения задачи с помощью соответствующих граничных условий. Для численного решения задачи (1) — (3) использовался метод динамической адаптации (H.A. Дарьин, В.И. Мажукин; 1987). В основу метода положен переход к произвольной нестационарной системе координат с переменными (<7,г), принадлежащей некоторому расчетному пространству fiit. В произвольной нестационарной системе координат неизвестными являются не только сеточные функции Т/, но и координаты узлов сетки л/. Для их определения используется уравнение обратного преобразования, представляющее собой дифференциальное уравнение в частных производных. Уравнение обратного преобразования составляется таким образом, чтобы скорость движения узлов расчетной сетки зависела от динамики решения уравнений, описывающих физические процессы. В произвольной нестационарной системе координат уравнение (1) принимает вид:

д(Ч>Н) ö(HQ) Ш

¡) т с) q д q

дгр Ь Q дт с)q

(4)

Ч> &У

(5)

д л" _ V i)q р

• Чя<Ч<Гк„, гьО, k-s,t,

где (5) - уравнение обратного преобразования, функция - <2 заранее неизвестна и подлежит определению. В качестве функции преобразования использовалась функция

<3=-Э —, обеспечивающая квазиравномерное распределение узлов в каждый момент ¿к/

времени, где О - коэффициент, определяемый через параметры задачи: геометрические

0.0- а)

Ч>.5 --1.0- 1 ' — Прямоугольный I аусспвсккй

-1.5-

-2.0-

-2.5- ч '

а)

Прямоугольный ;~ Га\*ссовский

Г>>

/

— Прямоугольный

I нусеонекий

\ /

. о)

' Прямоугольный Гауссовский

0,50.0- -

I, 10" с

Рис. 1. Скорости фронта испарения г>(,,(/), м/с Рис. 2. Толщина расплава Н„ 10"5м

—Трсугаиьный

ч вснрастаюший

\ _ _ Треугольный

ч ниспадающий

/ ч

/V ч

/ \

ч

-- V \

\

\

\

\

1,2 1,0 0,8 0.6 0.4 0,2 0.0 -0,2

■о .'Ч 2 4 6

1. 10"4 с

Я) 1\ -Треугольный

1 ( возрастающий

, } 1 реутопышй

$ , ниснадаюицм

у •

' 1 ;

-2 0 4 6

г. 10"1 с

а) л

■ -Треугольный

возрастающий

г . .. Треугольный

/ V нисиадаюаин:

И

.0)

1к — Треугольный вдарастакшшй .... Треугольный нислялаюошй

Рис. 3. Толщина расплава Н,, 10~5м

Рис.4. Скорости плавления 1А.Д?), м/с

размеры области и скорость движения границ О =-—-1, где Цп-длина области

V

в момент времени Г, - скорости движения фронта плавления и испарения

соответственно.

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

В качестве мишени рассматривались два металла с сильно отличающимися теплофизическими параметрами: медь и титан. Прямоугольная форма импульса выбиралась в качестве основополагающей. С ней сравнивались все остальные. При этом поглощенная мишенью интенсивность лазерного излучения составляла Сшг = 10" Вт/м2, а длительность -Г/.=210~4 с. Значения интенсивности 0(0 и длительности импульса выбирались такими, чтобы физически система успевала перейти в стационарное состояние как по температуре поверхности , так и по скоростям распространения

межфазных границ 1>1С<1'и,. Особенностью воздействия прямоугольных импульсов является наличие двух участков нестационарности, связанных с мгновенным изменением интенсивности - передний и задний фронты, соответствующие моментам включения и выключения импульса. Скачкообразное изменение интенсивности при включении приводит к формированию вблизи поверхности области с максимальными градиентами температуры, определяющими процесс плавления с максимальной скоростью на переднем фронте импульса и высокой скоростью затвердевания. Почти мгновенное прекращение испарения при медленно убывающем процессе плавления приводит к значительному, в 1,5 раза у меди и в 2 раза у титана, увеличению толщины жидкой фазы Я|Ш01э 85-Ю"6 м (Си), ЯЪ]1а1а 7,5-Ю"6 м (ТО (рис. 2). Сравнение толщин испаренного слоя Я, а 410-10"6 м и расплава // / //,е 5 для меди и Я, а 397-10"6 м , II % 111,3 53 для титана свидетельствует о том, что при использовании прямоугольных импульсов основная доля энергии импульса вкладывается в процесс испарения. Расчеты

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

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

Вторая глава посвящена разработке математической модели лазерного разрушения мишени в среде с противодавлением и построению вычислительного алгоритма решения задач радиационной газовой динамики (РГД), объединенной с гидродинамическим вариантом задачи Стефана. Постановка задачи осуществляется в рамках одномерной многофронтовой гидродинамической задачи Стефана, записанной для трёх фаз (твердая, жидкая, газообразная (пар, воздух)). Вычислительные особенности рассматриваемой задачи состоят в наличии многочисленных подвижных границ, среди которых межфазные и контактные границы, ударные волны и свободная граница. В основу вычислительного алгоритма положен конечноразностный метод динамической адаптации. Пространственное положение межфазных границ Г\(/), ударных волн Г„.,(/), Г,„.(г) контактной Г,,,,(г) и свободной Г,(г) границ изображено на рис. 5.

Твердая фаза Твердая фаза Жидкая фаза Пар Воздух Воздух Лазерное излучение <—-

■ С— I", Г,г I -.Ьо'к..' г, 1-> ■* А <-

Рис. 5. Пространственное положение межфазных границ

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

осуществляется следующей системой дифференциальных уравнений в частных производных (В. И. Мажукин, Г. А. Пестрякова; 1985):

д(> , Д(Н _0;

at tlx

а(ри) | | ар _0.

dt Эх дх

a(pe)td{p№) / Ои fdWr ^IW, tdG\

at вх ^ Их дх ax ax J

к i

w; = ffii/v'i/vi/,u. w, -iv; + и',-;

о -i

— + KL(p,T)G' = 0, G = G' +G~; ax

WT—l(T)y-. P = P(p,T), £ = e(p,T).

k~s,l,u,g.

Здесь p, и, e, T, P -плотность, газодинамическая скорость, внутренняя энергия, температура и давление материала соответственно; к и -коэффициент поглощения и спектральная плотность излучения плазмы; / <,,-плотность равновесного излучения; kl и G -коэффициент поглощения и плотность лазерного излучения;

lVr, IV, - тепловой и

радиационный потоки; G*,Wf- падающие на поверхность мишени лазерный и радиационные потоки; С,\ И]' - отраженные лазерный и радиационные потоки, А-коэффициент теплопроводности; и = Cos косинус угла вылета фотона. Индексы s, /, и, g обозначают принадлежность величин соответственно к твердой, жидкой, парообразной и воздушной фазам.

Для определения зарядового состава плазмы и коэффициента поглощения лазерного излучения использовалась нелинейная система уравнений Саха. Система (6) была дополнена начальными и граничными условиями, наиболее важными из которых являются граничные условия на подвижных границах, описывающие механизмы поверхностного плавления, испарения. Подвижные фазовые границы х - Г,, (/), х = Гк(/), контактная граница* = Г,,,(/) и фронт ударной волны дг = Г,,г(/) выделялись явным образом с помощью соответствующих граничных условий.

В новых переменных (д, т) система (6) примет вид:

а г а<?

(1т (V/ <)(/

д(уу^),а[р Ф'+б)] _ [рС'« , ¿ИУ,- , 'НУ, , а г а <7 ( а</ (1<? /

[I—~-+ук(1п>,р.Т)1,. = ;

ад

IV,* = Г Ги1.<1\'<1к, IV, = —^—;

а ■ V

а?

От

Р = Р(р,Т).е = е(р.Т), к=*.1.ъ8

где 6' = б" + (Г, V - якобиан обратного преобразования.

Функция преобразования для адаптации под большие градиенты решения определяется из принципа квазистационарности (В.И. Мажукин; 1993).

(у-1) я(р.г) а / |

а<?и-

1>-1р-к{дрдч ит дц а<г

ЛIV, | дО а<7 <iq

Я я

где ге - регуляризующая константа, ограничивающая снизу значение производной при её стремлении к нулю.

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

Система дифференциальных уравнений в частных производных аппроксимировалась разностными схемами, полученными интегро-интерполяционным методом (А.А.Самарский; 1977) на разнесенных сетках.

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

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

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

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

В качестве исследуемого материала рассматривалась алюминиевая мишень толщиной -30 мкм, находящаяся в воздухе с температурой 7V=273 К и давлением iV=105 Па. В качестве теплофизических параметров были использованы значения, соответствующие твердой и жидкой фазам алюминия. В качестве уравнения состояния для конденсированной фазы были использованы данные, приведенные в (K.S.Holian; 1986), а в газовой среде уравнения состояния Р = Р{р.Т), £ = с(/),Т) в табличном виде (Н.Н.Калиткин, И.В.Ритус, А.М.Миронов; 1983). Расчет коэффициентов поглощения производился с использованием модели Хартри - Фока - Слейтера.

Моделирование режимов лазерного воздействия производилось с длиной волны Я =0.8 мкм, с гауссовским по времени распределением интенсивности

0(/) = (/„ехр(-(//г)'), -оо</<ос длительностью rt=10"8-40"15 с и пиковой интенсивностью Go = 1012- Ю20 Вт/м2. В расчетах плотность энергии J = С0г, менялась в диапазоне Ю4 £,/ £ 5-Ю5 Дж/м2, поглощательная способность поверхности А =10%. Результаты математического моделирования показали, что плазма в испаренном веществе образуется для J >3 105 Дж/м2 (пороговое значение для плазмы Jmax~(3^5) Ю5 Дж/м'). Так как плазма замедляет процесс испарения, основное внимание уделялось доплазменным режимам.

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

и, как известно, при г-»0 скорость фронта плавления 5 ~ Г1" —> ос. Однако

реальная скорость плавления всегда остается конечной величиной. Одним из эффективных механизмов ограничения скорости 1'!( является гидродинамика расплава. Её влияние проявляется, в частности, в зависимости температуры плавления Тт = Тт (/') от давления Р, на плавящейся поверхности, величина которого определяется значением скорости 1\,. В итоге рост скорости I',, приводит к росту давления Р, и температуры на межфазной границе Т,,= Т„ (Р,), а увеличение температуры плавления Т„~Тт{Р1) ограничивает рост скорости . Типовой пример пространственной структуры решения при ультракоротком г = 10 фс , / = 104 Дж/м2 воздействии на момент формирования ударных волн в твёрдой фазе и газовой среде Г = 80 фс представлен на (рис.6).

Объёмный механизм поглощения излучения и поток вещества через межфазную границу х = Г, формируют в твёрдой фазе приповерхностный температурный

максимум.

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

Так как между максимальными значениями скоростей и ии, выдерживается

соотношение к, »и,,., то значение максимума температуры в твердом теле намного превосходит максимум температуры в жидкости.

X, 10-6.4

Рис. 6. Пространственный профиль температуры в момент I =78.5 фс для режима воздействия т = Ю фс, Сп = 10" Вт/м2

Основная особенность пороговых для абляции режимов г, < 10 '"с и

10' < ,/ < 104 Дж/м2 состоит в том, что толщина жидкой фазы не превышает 100 + 200 А и она прозрачна для лазерного излучения. Энергия лазерного импульса почти полностью поглощается твердой фазой, вызывая её значительный перегрев (Т1тк ~ 2000ЙО-Запасенная энергия в основном затрачивается на нагрев и плавление твердой фазы. Поэтому процесс плавления протекает более эффективно, а роль испарения при низких плотностях энергии 10' < .7 < 104 Дж/м2 незначительна. Для усиления процесса абляции в пикосекупдном диапазоне необходимо увеличение поглощаемой плотности энергии до /> Ю4 Дж/м2.

В диапазоне ультракоротких импульсов (И) <т < ¡0 ) секунды при 104 й./5105 Дж/м2 быстрый нагрев мишени, высокая скорость плавления и относительно медленная передача энергии за счет теплопроводности приводят к существенно большему перегреву твердой фазы по сравнению с длинными наносекундными импульсами. Максимальная скорость фронта плавления при г, =10"|3*10"и с достигает скорости звука в твёрдом теле, которая для алюминия составляет ~ 6.26 103 м/с.

Приближение скорости к скорости звука, с одной стороны, создаёт условия для возникновения ударной волны в твёрдой фазе, а с другой - способствует её максимальному перегреву, (см. рис. 6.). На момент образования / =80 фс ударных волн в газовой и твёрдой фазах, максимум температуры в твердом теле при т = 10 фс и

У = 104 ДжЛг, находится на глубине -200 Л от границы Г,. (/) и составляет 6-10' К. Степень перегрева Т!т>у.1Тм достигает величины ~ 6.5. Соответственно скорость, температура и давление на границе раздела фаз Г5.,(/) таковы: =4■ 10! м/с, Т »Тт1) = 1.5-Ю3 К, и Р5 = 5 ■ 10ю Па. При полученных значениях давления возбуждается сильная ударная волна в твердой фазе.

Рис. 7. Пространственный профиль температуры в момент ! -2,5 фс для режима воздействия г = 1 пс, С0 = 10" Вт/м"

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

примерно одинаковым и равняется / = (2.5 + 3) не. Приповерхностные максимумы температуры полностью исчезают к моменту Г ~ (2.5 + 3) не (рис. 7).

Для наносекундных импульсов плазма в испаренном веществе не образуется, если плотность излучения J < 3 105 Дж/м2. При тех же значениях У для фемтосекундных импульсов термическая плазма возникает после окончания воздействия.

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

Основные результаты диссертационной работы:

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

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

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

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

Установлены пороги плазмообразования в испаренном веществе и окружающей газовой среде. Для наносекундных импульсов плазма образуется в испаренном веществе для J >3 105 Дж/м2. Для фемтосекундных импульсов плазма возникает в воздухе после окончания лазерного воздействия и носит термический характер.

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

1. Mazhukin V.l., Lobok M.G., Chichkov B.N. Modeling of fast phase transitions dynamics in metal target irradiated by pico and femto second pulsed laser // Applied Surface Science. 2009. № 255. P. 5112-5115.

2. Мажукин В.И., Мажукин A.B., Лобок М.Г. Математическое моделирование динамики фазовых переходов и перегретых метастабильных со стояний при нано -фемтосекундном лазерном воздействии на металлические мншени // Математическое моделирование. 2009. Т. 21, № 11. С. 99-112.

3. Mazhukin V.l., Mazhukin A.V., Lobok M.G. Comparison of nano- and femtosecond laser ablation of aluminium // Laser Physics. 2009. V. 19, № 5. P. 1169 -1178.

4. Mazhukin V.l., Lobok M.G., Smurov I. Transient effects in pulsed laser irradiation // Applied Surface Science. 2007. V. 253. P. 7744 - 7748.

5. Лобок М.Г., Мажукин В.И. Влияние временного профиля импульсов на процессы лазерного воздействия // Математическое моделирование. 2007. Т.19, №9. С. 54-78.

6. Lobok M.G., Chuiko М.М. Dynamic adaptation method for the solution of Stefan problem // Proceedings of the III International Conference on Adaptive Modeling and Simulation. Göteborg, 2007. P. 150-153

7. Лобок М.Г. Математическое моделирование влияния теплофизических параметров и временной формы импульса на процессы лазерного плавления и испарения металлов // Ломоносов-2004: Труды Международной научной конференции, М., 2003. С. 283-284.

Подписано в печать: 21.10.2013 Объем: 1,1 усл.п. л. Тираж: 100 экз. Заказ № 930 Отпечатано в типографии «Реглет» 101000, г. Москва, Пл. Мясницкие Ворота д.1, стр.3 (495) 971-22-77 www.reglet.ru

ИНСТИТУТ ПРИКЛАДНОЙ МАТЕМАТИКИ ИМ. М.В. КЕЛДЫША РАН

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

ЛОБОК Максим Геннадьевич

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ИМПУЛЬСНОГО ЛАЗЕРНОГО ВОЗДЕЙСТВИЯ НА МЕТАЛЛЫ

Специальность 05.13.18 - математическое моделирование, численные методы и комплексы программ

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель: д.ф.-м.н., профессор Мажукин В.И.

Москва - 2013

Оглавление Стр. Введение .......................................................................................................................4

1. Режимы лазерного воздействия...........................................................................4

2. Проблема Стефана. Основные математические модели...................................5

3. Динамическая адаптация...................................................................................11

4. Лазерное воздействие.........................................................................................14

Глава 1. Влияние временной формы лазерного импульса больших длительностей и теплофизических характеристик среды на динамику фазовых переходов...................................................................................................................21

1.1. Введение...........................................................................................................21

1.2. Постановка задачи...........................................................................................23

1.3. Теплофизические и оптические характеристики и параметры...................25

1.4. Метод и алгоритм численного решения........................................................25

1.5. Анализ результатов моделирования..............................................................32

1.6. Выводы.............................................................................................................54

ГЛАВА 2. Вычислительный алгоритм для задач лазерной абляции в режимах с возникновением плазмы.....................................................................................56

2.1. Введение...........................................................................................................56

2.2. Постановка задачи...........................................................................................57

2.3. Произвольная нестационарная система координат......................................67

2.4. Выбор функции адаптации.............................................................................70

2.5. Разностная аппроксимация.............................................................................73

2.6. Выводы.............................................................................................................78

3 Стр. Глава 3. Математическое моделирование динамики фазовых переходов при нано-фемтосекундном лазерном воздействии на металлические мишени .....................................................................................................................79

3.1. Введение...........................................................................................................79

3.2. Постановка задачи...........................................................................................82

3.3. Вычислительный алгоритм.............................................................................83

3.4. Входные данные...............................................................................................83

3.5. Результаты моделирования.............................................................................85

3.6.Вывод ы..............................................................................................................99

Выводы ...................................................................................................................100

Введение

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

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

1. Режимы лазерного воздействия

Большинство технологических приложений импульсной лазерной обработки материалов связано с началом фазовых превращений 1-го рода: плавления-затвердевания и испарения-конденсации. Соотношение глубин расплава и испарённого слоя играет важную роль, определяя качество таких технологических операций как лазерное сверление, резка и сварка. В этих операциях размеры расплавленной зоны и качество отверстия зависят, как от режима воздействия, определяемого длиной волны, длительностью импульса и пространственно-временным распределением интенсивности в импульсе, так и от теплофизических и оптических свойств обрабатываемого материала [1]. При разработке конкретного технологического приложения влияние каждого из факторов должно быть исследовано и охарактеризовано. Определение их роли связанно с моделированием пространственно-временных распределений температурных полей и динамики фазовых переходов. Для этих целей использовался ряд математических моделей [2] - [4], часть которых

допускает аналитические решения [5], [6], но подавляющее большинство решается численно [7].

2. Проблема Стефана. Основные математические модели

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

2.1. Классический вариант задачи Стефана. Плавление-кристаллизация

Классический вариант задачи Стефана, сформулированный для переходов типа плавления или затвердевания, сводится к уравнению теплопроводности в области с заранее неизвестной подвижной границей которая разделяет твердую ¿■и жидкую С фазы.

На межфазной границе выполняется так называемое дифференциальное условие

Стефана, которое имеет простой физический смысл: скорость движения фазовой границы и!( (/) определяется разностью потоков тепла, и пропорциональной поглощаемой или выделяемой на этой границе объемной теплоты фазового перехода рхЬт.

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

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

Между тем в классическом варианте задачи Стефана неравновесность фазового превращения явно не фигурирует. Подобный подход и его различные модификации широко используются для списания фазовых переходов при импульсном воздействии концентрированных потоков энергии [8,9] и имеют, безусловно, определенную область применимости, однако вопрос о ее границах для ряда практически важных случаев остается пока открытым. Использование предположения о постоянстве температуры фазовой границы оправдано, например, для переходов твердое тело-жидкость при относительно малых ин-тенсивностях излучения, которым соответствуют медленные изменения температуры и малые значения скорости (/).

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

2.2. Энтальпийная постановка

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

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

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

23. Неравновесная модель

Использование предположения о постоянстве температуры фазовой границы может быть оправдано, например, для переходов твердое тело - жидкость при относительно

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

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

2.4. Однофазный вариант задачи Стефана. Испарение

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

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

Для описания процесса поверхностного испарения обычно применяется однофазный вариант задачи Стефана [15], в котором температура плоского фронта перехода уже не является постоянной, а слабо, логарифмически, зависит от скорости фронта vCv(t). Уравнение теплопроводности в задаче Стефана дополняется тремя законами сохранения: массы, импульса и энергии и двумя дополнительными соотношениями на внешней стороне кнудсеновского слоя. Дополнительные соотношения характеризуют кинетику фазового перехода и степень неравновесности процесса испарения. Дополнительные соотношения обычно получают при некоторых предположениях о неравновесной функции распределения внутри кнудсеновского слоя. Имеется несколько способов аппроксимации кнудсеновского слоя [16-29].

2.5. Гидродинамический вариант задачи Стефана

Одномерное гидродинамическое движение в конденсированных средах может возникать в результате изменения плотности вещества под влиянием теплового расширения в области прогрева и изменения объема вещества при фазовых превращениях.

Так при сильном перегреве скорость фронта в классическом варишгге задачи Стефана может формально обращаться в бесконечность. Это следует, например, из точного решения одномерной задачи [30] о плавлении твердого полуограниченного образца,

начальная температура которого Т5 > Тт превышает температуру плавления. В этом случае в начальный момент времени скорость г^Д/) неограниченно велика, а при больших перегревах это решение вообще теряет смысл. Корректное математическое моделирование таких физических ситуаций требует привлечения полной системы гидродинамических уравнений [31], в то время как в классическом варианте задачи Стефана рассматривается лишь уравнение переноса энергии, без учета изменения плотности и давления.

Гидродинамический вариант задачи Стефана в одномерном приближении используется также для описания проблем, связанных с генерацией и распространением оптоаку-стических сигналов [32]. Гидродинамический вариант представляет собой краевую задачу для полной системы уравнений гидродинамики, дополненной уравнением энергии в области, разделенной подвижной фазовой границей Г,((/).

Данная система дополняется соответствующими уравнениями состояния и граничными условиями. На фазной границе выписываются соотношения характерные для одного из рассмотренных вариантов задачи Стефана.

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

3. Динамическая адаптация

Проблеме построения расчетных сеток при решении задач математической физики в последние годы уделяется большое внимание [33-95]. Существует мнение [36,44], что окончательное разрешение многих проблем численного решения задач математической физики следует ожидать не только от улучшения способов разностной аппроксимации уравнений в частных производных и усовершенствования алгоритмов решения сеточных уравнений, так и от правильного выбора расчетной сетки. Точность решения уравнений в частных производных зависит от, согласованности распределение узлов сетки с особенностями искомого решения. Из двух решений одной задачи, полученных на двух различных сетках с одинаковым числом фиксированных узлов, меньшая погрешность будет достигаться на сетке с более оптимальным распределением их по отношению к искомому решению. Принцип оптимального распределения узлов положен в основу методов построения адаптирующихся к решению сеток. По признаку имеющихся особенностей решения наиболее важные задачи, требующие применения адаптивных сеток, удобно разбить на следующие классы [36,40]:

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

• задачи с большими градиентами вблизи границы области: к ним относятся задачи типа по�