Численное моделирование нагрева металлической мишени фемтосекундными лазерными импульсами в присутствии газа Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 13, № 3, 2008

Численное моделирование нагрева металлической мишени

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

*

в присутствии газа

Н. М. Булгакова Институт теплофизики СО РАН, Новосибирск, Россия

В. П. Жуков

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

We present the results of numerical modeling of a thermal exchange between a metal target (platinum) and an ambient gas (argon) under the action of the femtosecond laser pulses on the target. Dynamics of gas motion is described by the Navier—Stokes equations and the heating process is governed by the heat flow equation which takes the phase transitions into account. It is shown that gas breakdown in the focus of the laser beam leading to the subsequent heat transfer from the hot gas to the target can explain the experimentally observed effect of an abnormal heating of the metal targets. A finite-difference algorithm for solving the problem is described and we present the analysis of the physical processes which take place under the interaction of the ultra short laser pulses with metals in the presence of the ambient gas.

Введение

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 06-08-01196 и № 08-01-00264-а).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

1. Анализ физических процессов

В экспериментах [5-7] металлическая мишень, окруженная газом (варьируются химический состав и давление газа), облучалась лазерными импульсами длительностью 65 фс, энергией до 1.5 мДж и длиной волны 800 нм. С помощью калоритмической техники измерялось тепло, оставшееся в мишени после облучения. Экспериментальные данные показывают, что в вакууме и газе при малых энергиях лазерного импульса нагрев мишени соответствует величине, определяемой коэффициентом отражения металла Rf, как Ег = (1 — Rf)Е¿, где Ег — остаточное тепло в мишени, а Е^ — энергия лазерного импульса. В случае абляции в фоновом газе при увеличении мощности излучения остаточное тепло в мишени может увеличиваться более чем в 2 раза по сравнению с величиной Ег. Варьирование состава газа (использовался аргон, неон, гелий при давлении до 1.08 атм и комнатной температуре) убедительно свидетельствует о том, что эффект аномального нагрева мишени зависит от потенциала ионизации газа. Чем выше потенциал ионизации молекул газа, тем большая энергия лазерного излучения необходима для возникновения этого эффекта.

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

Причиной ионизации газа в нашем случае может служить только многофотонная ионизация. Изменение числа электронов со временем в этом случае описывается формулой [15]

дпе/ог = АН1к (г),

где N — плотность нейтральных атомов, А и к — константа и порядок многофотонной ионизации соответственно, 10 — мощность лазерного излучения. Оценки показывают [16], что при параметрах экспериментов [5-7] мощность излучения примерно 3 Дж/см2 достаточна для 100 %-й однократной ионизации аргона в области фокуса лазерного луча.

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

2. Математическая постановка задачи

Задача об энергообмене нагретого излучением аргона и металлической мишени моделируется следующим образом. Течение газа описывается уравнениями газодинамики с учетом теплопроводности и вязкости, а нагрев металлической мишени (платина) — уравнением теплопроводности с заданными коэффициентами теплопроводности и теплоемкости. Эти коэффициенты зависят от фазового состояния мишени. Испарением мишени в данной работе пренебрегается, а металл, нагретый выше точки плавления, считается неподвижным. Основные уравнения имеют вид

/ д v \

P\~dt +(vv)vj = —vp + div П, (1)

dp + div (pv) = 0, (2)

dU dV

— + div (vU) = —p div v + div (xVT) + (7r)ra>mдхт, (3)

^ (dVn dVm 2 (n)n,m = V i д--+ д--3"ramdiv v

\dxm dxn 3

Здесь v, p, p и v — скорость, плотность, давление и вязкость газа; U и T — внутренняя энергия и температура газа соответственно.

Задача решается в цилиндрической (r—z) геометрии. Предполагается, что газ занимает область 0 < r < R, — Zg < z < 0, а платина — область 0 < r < R, 0 <z< ZPt. Соответственно, уравнения (1) и (2) решаются в области — Zg < z < 0, а уравнение для энергии (3) — в области — Zg < z < Zpt. В области, занимаемой платиной, скорость в (3) полагалась равной 0. Связь между давлением, температурой и внутренней энергией газа соответствует идеальному газу с показателем адиабаты y:

p =(Y — 1)U;

p = po(p/po)(1+ T/T- c). (4)

Здесь p0 = 1.013 • 105 Па — давление, соответствующее 1 атм; p0 — плотность газа при 1 атм и 0 °C. Температура измеряется в градусах Цельсия, T-с = 273.15 K.

В области, занятой платиной, решается уравнение (3), в котором полагается v = 0. Если температура не равна температуре фазового перехода, имеем соотношение U = J C(T)dT + const, где C — объемная теплоемкость. В расчетах теплоемкость полагается постоянной в пределах одной фазы. Если температура равна температуре фазового перехода, T остается постоянной при изменении внутренней энергии на величину, равную скрытой теплоте перехода (задача Стефана).

Заметим, что при переходе через границу z = 0 энергия может изменяться скачком, поскольку она может отсчитываться от разных уровней в газе и мишени, однако температура и тепловой поток —xdT/dz должны быть непрерывными. Это условие является граничным условием для уравнения (3) при z = 0. Остальные граничные условия таковы. При z = 0 компоненты скорости газа равны 0 (условие прилипания). На периферийных границах в области газа ставятся условия проскальзывания:

Vz = 0; dVr/dz = 0 при z = — Zg,

К = 0; дУг/дг = 0 при г = R,

а температура на этих границах полагается фиксированной и равной фоновой температуре системы Т0 = 20 °С. На периферийных границах в области, занятой платиной (г > 0), ставятся условия отсутствия потока тепла (нормальная производная температуры равняется 0). Заметим, что периферийные границы отодвигаются далеко от области нагрева, и конкретный вид граничных условий на них не играет роли. Граничные условия на оси имеют обычный вид

К = дУх/дг = ди/дг = др/дг = дТ/дг = 0.

Начальные условия соответствуют неподвижному (v = 0) газу при постоянной плотности. Газ в области фокуса лазерного излучения размером RL по г и по г предполагается мгновенно нагретым. Распределение температуры имеет вид

Т(г = 0) = То + Ть ехр(—г2/R2L — г2/Х\), где TL определяется из выражения

^2п / ехр(—г2^ — г2/Х1 )ЫЫг = ^^^ + Р^^(5)

7 — 1 ро Т°с ] 7 — 1 ро Тс

Здесь интеграл берется по всей области, занимаемой газом; Т0 = 20 °С — фоновая (комнатная) температура; Тюп = 15.76 • 11604(7 — 1) К — температура, соответствующая энергии ионизации, равной в случае аргона 15.76 эВ (1 эВ соответствует 11 604 К); Rf — коэффициент отражения платины; QL — энергия лазерного импульса на единицу площади пятна облучения на поверхности мишени; пR2LQL — полная энергия лазерного импульса. Выражение (5) соответствует тому, что выделившейся в газе энергии лазерного импульса достаточно для ионизации аргона в объеме 'кRLZL и, кроме того, в образовавшейся плазме поглощается дополнительно некоторая доля (Ргс) от энергии лазерного излучения, отразившегося от платиновой мишени Rf пRLQL за счет эффекта интерференции падающего и отраженного от мишени излучения [16]. Величина Ргс является параметром задачи.

Распределение внутренней энергии в платине в начальный момент времени имеет вид

и = ир ехр(—^г — г2 /'К^).

Здесь ^ — коэффициент поглощения излучения. Величина и0р4 вычисляется из условия равенства полной энергии, заключенной в платине и0Рь2п [ [ ехр(—^г— г2/^^)гйгйг (интегрирование проводится по всему объему, занимаемому платиной), энергии, поглощенной мишенью (1 — Rf)пRLQL.

3. Значения параметров

Задача решалась для платиновой пластины толщиной Хр^ = 1 мкм = 10-6 м. Мишень облучалась лазерным импульсом с интенсивностью 3 Дж/см2 и размером пятна облучения RL = 50 мкм. Толщина зоны пробоя в фокусе лазерного луча перед мишенью составляла ZL = 100 мкм. Размер расчетной области по радиусу и размер области, занимаемой газом, по г изменялись от 220 до 3000 мкм по мере того, как движение газа охватывало все большее пространство.

Начальная температура системы платина—аргон составляла 20 °С. Аргон находился под давлением 1.08 атм. Показатель адиабаты аргона равен 7 = 5/3. Зависимости теплопроводности и вязкости аргона описывались выражениями хлг = 0.016(1 + Т/ТоС)1/2 Вт/(м-К) и V = 0.00221(1 + Т/ТоС}1/2 кг/(м-с).

Теплофизические и оптические свойства платины зависят от температуры, и их значения различаются в различных справочных изданиях. В настоящей работе взяты средние, постоянные для различных фаз значения: коэффициент поглощения излучения ^ = 50 мкм-1, коэффициент отражения Rf = 0.65, теплопроводность твердой и жидкой фаз составляет 74 и 55 Вт/(м-К) соответственно, теплоемкость твердой и жидкой фаз — 2.992-106 и 4-106 Дж/(м3-К) соответственно, температура плавления платины равна 1769.3 °С, теплота плавления — 2.37-109 Дж/м3, температура кипения платины — 3827°С, теплота испарения — 4.9-1010 Дж/м3.

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

4. Нормировка

При проведении расчетов в качестве единиц измерения используются следующие величины:

— длина — z0 = 1 мкм = 10-6 м,

— время — t0 = 1 нс = 10-9 с,

— скорость — z0/t0,

— плотность — р0 = 1.7839 кг/м3, что соответствует плотности аргона при T = 0 °C и давлении 1.013 • 105 Па (1 атм),

— давление и внутренняя энергия — p* = p0z0/t2. Температура измеряется в градусах Цельсия.

Различные коэффициенты в безразмерном виде связаны с размерными коэффициентами (знак dim) следующим образом:

— теплопроводность х = Xdimt(]/(P0z0) = XdimW(p*z0),

— вязкость v = vdimt0/(p0z2),

— объемная теплоемкость C = Cdim/p*,

— коэффициент поглощения ^ = ^dimz0,

— плотность лазерного излучения QL = QLim104/(p*z0) (коэффициент 104 связан с тем, что интенсивность Ql задана в Дж/см2, а не в Дж/м2).

В такой нормировке вид всех уравнений и соотношений, за исключением (4) и (5), сохранится. Формулы (4) и (5) примут вид

p = вр(1 + T/T-C), (6)

-^т-^2п /exp(-r2/R| - z2/ZL)rdrdz = ^ь + PrcRfnR£Ql,

Y - 1 P0 J°c J Y - 1 P0 T°C

где в = p0/p*. Заметим, что значение давления, измеряемого в атмосферах, и безразмерное давление связаны между собой соотношением patm = p/в.

В цилиндрической геометрии уравнения (1)-(3) принимают следующий вид:

'ВУг ЛГ6Уг ЛгдУг\

+ + У.- г *

дь

дг

дг

др 1 д(гп„) дп.г п^

--о--1---о--г

дг г дг

дг

дУ дУ* т 6У2

Р\— + Уг— + Уг

дь

дг

дг

др дп.

дг дг

^ + 1 д (гПгг)

г дг

пг

V --

2 дУг + 4 дУ 2 Уг

3 дг 3 дг 3 г

п

Пг

V — -

2 дУг 2 дУг + 4 Уг

3 дг 3 дг 3 г 4 дУг 2 дУ. 2 Уг

V -

3 дг 3 дг 3 г

Пгг = V

дУг дУг

г + г

дг дг

др + а1у (рУ) = о,

ди

— + (у и) = -р у + (хут ) + дь

дУг дг

Уг

дУг

Пг

+ --+ пгг—--1--

дг

V

а1у { = 1 дМ) + д£, уТ = (дТ,дТ

г дг дг дг дг

(7)

(8) (9)

10)

11) 12)

13)

14)

15)

5. Конечно-разностная схема

Для численного решения поставленной задачи использовалась конечно-разностная схема расщепления типа схемы стабилизирующей поправки [17] на сдвинутой адаптивной квазиравномерной [18] регулярной сетке. Давление, плотность, температура, внутренняя энергия, коэффициенты вязкости, теплопроводности и теплоемкости и компоненты тензора вязких напряжений пгг, , пгг вычислялись в центре ячеек разностной сетки (г, ] соответствуют точке с координатой гг, г^), скорости — на гранях (УГг+1/2,^, У ^+1/2), а компонента пгг — в узлах г + 1/2, ] + 1/2.

Узлы сетки располагались таким образом, что гг=0 = 0 (ось цилиндра), г^=^о+1/2 = 0 (граница аргон—платина). Вводились также фиктивные узлы с координатой г-1/2 = —г1/2. Это позволяло простым и естественным образом аппроксимировать изучаемые уравнения. Уравнения (8), (13), (14) решались в том числе и на оси цилиндра, т. е. точки гг=0 не являлись границей. Узлы г-1/2 были граничными для Уг, при этом Уг -1/2 = Уг 1/2.

Пространственные производные в уравнениях (7)-(15) аппроксимировались центральными разностями со вторым порядком обычным образом. Например:

Уг

дУг

дг

г+1/2,3

(Уг г+1,^+1/2 + Угг+1,]-1/2 + г,^+1/2 + У г,.7-1/2)( У г+1/2 ,.?+1 — Уг+1/2,.7-1)

4(^ + 1 — ^-О

г

г

д (ЦУ)

дг

г;

+ — + Цг—1; ^г; —1/2

2(г.?+1/2 - ^^-1/2)

Отметим некоторые особенности аппроксимации членов вида г 1д(г/^Т)/дг и ^/г в точках гг, связанные с наличием оси г = 0. Для них имеем

1 д(г/Уг) \ г дг /

1 гг+1/2 (/г+1 + /г)Кг+1/2 — гг—1/2(/г + /"г—1) К г—1/2

г=0

1 д(г/Уг) \ г дг у

2(г,+1/2 - гг-1/2) (/0 + ЛЖ 1/2

г=0

к

г

гг+1/2;

г1/2 ) + ^г г—

(16)

1/2,;

¿=0; Уг

г

2гг

К

¿=0;

г 1/2,;?

г1/2 '

Еще одной особенностью является вычисление производной дУг/дг вблизи границы аргон—платина г = 0. В уравнении (7) необходимо вычислять эту производную в точке ] = ^0 — 1/2, а в выражении (12) — в точке ] = В этих случаях функция УТ(г) аппроксимируется полиномом второй степени по г, коэффициенты которого находятся из условий К-(г = 0) = 0 (граничное условие), К1/2) = 1/2, К(2^—3/2) = 3/2. Производная д % /дг полагается равной производной от этого полинома в соответствующих точках.

Заметим, что для вычисления необходимо знать вязкость V на границе аргон— платина. Поскольку V является функцией температуры, то возникает вопрос о значении Т на этой границе. Этот вопрос тесно связан с вопросом об аппроксимации выражения ^у (ХУТ) в случае разрывных коэффициентов теплопроводности и резкого изменения шагов пространственной сетки, которое может иметь место на границе аргон—платина при г = 0. Подчеркнем, что уравнение (14) решается сквозным образом, как в области, занимаемой газом, так и в области, занимаемой платиной. Отличия состоят в том, что в последнем случае полагается v = 0, а также изменяется связь между Ц и Т. Заметим также, что коэффициент х терпит разрыв на границе различных фаз в платине.

Для построения разностного аналога ^у (хVT) запишем приближенное равенство на границе двух ячеек г;+1/2 = 0

Т;

/дТ\

;+1 = Т;+1/2 + ( ) (г?+1 — г;+1/2), V дг / ;+1/2+0

(17)

Т; = Т;+1/2 — i —

дг) (г;+1/2 — )-дг / ;+1/2_0

(18)

Здесь (дТ/ дг);+1/2±0 — правые и левые производные на границе двух ячеек. Из непрерывности теплового потока имеем

Т

Х;+1/2+01 ^

;+1/2+0

х ;+1/2—0

дТ

дг

;+1/2—0

г

г

Полагая постоянство х на промежутках zj-i/2 < z < zj+1/2, имеем

дТ\ (дТ\

М Tz) = j dZ • (19)

j+1/2+0 \dZ/ j+1/2-0

Из (17)-(19) для потока тепла = — xdT/dz в узлах zj+1/2 имеем

_ Xi.j+iXi.j

+1/2 = — х (Z Г"

Xij+1 VZi,j+1/2 —

Соответственно

^j+v2 =— х (z ; ^TtxTZ z тт(Ti,j+i — ^

Х i,j+1 VZij+1/2 — zj ) + Xij (zj+1 — zj+1/2)

д ( <9T^ = qzi,j+1/2 — 9zi,j-1/2 X

i,j Zj+1/2 — Zj-1/2

Аналогичные формулы имеем для потока qri+1/2,j. Для r-1Sqr/5г при r = 0 имеем выражения, аналогичные (16).

Из (17)-(19) следуют выражения для температуры на гранях ячеек, в частности, на границе аргон—платина:

T = Х Tj (zj+1 — zj+1 / 2j + Xj+1Tj+1 (Zj+1/2 — zj ) Xj+1 (Zj+1/2 — zj j + Xj (zj+1 — zj+1/2)

Переход со временного слоя tn на слой ¿""+1 = tn + т осуществляется следующим образом.

1. Зная значения pn, Un, vn, вычисляются Tn, pn, vn, Пn, xn, Cn.

2. Решается уравнение для плотности

P(1) — pn

P-^ = L1pn + L2pn,

(tL1 — 1)(р(2) — pn) = —(р(1) — pn), J2 — 1)(pn+1 — Pn) = — (P(2) — P"

= 1 д (rVnp) = d(Vnp)

L1P =----j l2P =---•

r dr dz

3. Далее вычисляется внутренняя энергия

U(1) — Un

■ = L1Un + L2Un — pndiv vn + Q

т

№ — 1)(и(2) — ип) = —(и(1) — ип), (тЬ2 — 1)(ига+1 — ип) = —(и(2) — ип).

В области, занятой аргоном, выражения для и Ь2 имеют следующий вид:

(7 — 1)Т°С 1 д ( пд(и/рга+1Л 1 д(гУгпи) 2 1 д/ г ди

1 в r 5r у X dr j r dr + (^ ) r Sr \vpn+1 5r

т U (7 — 1)Tqc д (д(U/Pn+1)\ d(VnU) + (y 1)2Un д ( z 5U" L2U =-в-dz Iх dz ) — + T(7 — 1) U д^ Vpn+T ^

Последний член в выражениях для Ь и ¿2 добавлен для лучшей численной устойчивости расчета по отношению к возмущениям звукового типа [19]. В области, занятой платиной, имеем

ъи = ~ (тх'^ 1 , ¿2и = # О9(и'С-)

г дг \ 5т у дг \ дг

Здесь используется теплоемкость платины в твердом состоянии С8. Это связано с тем, что в изучаемом течении наибольшее время платина находится в твердом состоянии, а значит, конечно-разностная схема должна обладать наилучшей устойчивостью именно в этом случае. Заметим также, что в ячейках сетки и ^ + 1, примыкающих к границе аргон—платина, оператор ¿2 должен быть записан с учетом того, что в его построение вовлекаются как ячейки, содержащие газ, так и ячейки, содержащие платину. Также должен быть изменен последний член в выражении для ¿2 в ячейке (аргон) с учетом граничного условия для скорости ^^/2 = 0. Полученные выражения здесь не приводятся в силу их громоздкости.

4. На следующем этапе решается уравнение для V:

уг(1) - V' 1 <%'+1 /

р'+1 дг р'+1:

№ - 1)(К(2) - V') = -(К(1) - V'), (т^2 - 1)(у+1 - V') = -(К(2) - V'),

у = 4^^ ( сЖ\ - 4^К - уядК 1 г 3 р'+1 г дг V дг У 3 р'+1 г2 г дг

ы = и ы1 - V" №

р'+1 дг V * дг '

Здесь под /г подразумеваются члены уравнения (7), связанные с тензором вязких напряжений, которые не вошли в операторы ¿1, ¿2. 5. Далее решается уравнение для V*:

^(1) - V' г т/п _ т/я 1 /'

* - ¿1^' + ¿2^' - „ч о + *

* ' ' * р'+1 ^ р'+1'

(тЬ - 1)(К(2) - к') = -(к(1) - V'), № - 1)(К'+1 - V*') = -(К(2) - V*'),

^ = Л- (г*' ^ 1 - ^'+1 ,

р'+1 г дг \ дг у г дг

Ь у = (^ 1 - у'

2 * 3 р'+1 дД * дг '

Под /* подразумеваются члены уравнения (8), связанные с тензором вязких напряжений, которые не вошли в операторы Ь1, Ь2.

Остановимся на используемой сетке. Размер области, занимаемой платиной по г, не менялся при проведении расчета и составлял = 1000 мкм. Размеры расчетной области по радиусу Я и размер области, занимаемой газом по оси г (величина увеличивались по мере распространения волны в газе от 300 до 3000 мкм по обеим координатам.

Сравнение результатов расчетов при разных К и Zg показывает, что увеличение размера расчетной области до значений более 3000 мкм не приводит к существенным изменениям. Используется квазиравномерная [18] сетка гг, г^-. Расчеты (см. ниже) показывают, что в газе наибольшие градиенты параметров имеют место в области расходящейся сферической волны и узкого по г слоя в окрестности границы аргон—платина г = 0. Наличие этого слоя связано с большой разностью температур мишени и газа. Кроме того, большие градиенты различных функций по г возникают и существуют длительное время в области, занимаемой платиной в окрестности г = 0, — это связано с тем, что энергия излучения поглощается в поверхностном слое платины. В областях больших градиентов шаг сетки устанавливался меньше, чем в остальной области. В процессе расчетов сетка адаптировалась таким образом, чтобы обеспечить необходимое количество точек в областях больших градиентов различных функций. Для представленных ниже результатов сетка менялась в процессе расчета около 10 раз. Для перехода с одной сетки на другую была создана программа пересчета значений различных функций.

Сетка по г в расчетах строилась из следующих соображений. Выделялись три области: область невозмущенного газа г2 < г < К, область фронта ударной волны г1 < г < г2 и область за фронтом волны г < г1. Внутри каждой области шаг сетки был близок к постоянному. Переход от одного шага к другому осуществлялся плавно на расстоянии порядка 8. Если 8 ^ 0, то шаг сетки менялся резко при переходе от одной области к другой. Чем больше 8, тем изменение шага плавнее.

Конкретно в расчетах при построении сетки по г использовались следующие формулы:

Лг С С С

^ = ,()+/( ), /(г) = 1 + С1 Ш((г — п)/81) + С2 Ш((г — г2)/82) (20) ^ / (г) + / ( —г) 81 82

или

F(ri) = iC, F(r) = 2r + C In fCh(; - + C2ln (Ch(r - - (21)

Vch(r + ri)/^lV \ch(r + ^/¿2)/

Значение C находится из условия ri0 = R, т. е. C = F(R)/i0. Уравнение (21) решается методом хорд. Формула (20) соответствует тому, что при ^ 0 и i0 ^ то шаг

сетки h = dr/di на промежутках r < r1, r1 < r < r2 и r2 < r < R постоянен и пропорционален числам h0, h1 и h2 соответственно. Для этого необходимо положить C1 = (h2/h1 — h2/h0)i1, C2 = (1 — h2/h1)i2. Для того чтобы при приближении к оси сетка стремилась к равномерной, в (20) для dr/di взята четная по r функция. Для полуцелых точек имеем уравнение F(ri+1/2) = (i + 1/2)C. Аналогично строится сетка по z.

В течение нескольких первых микросекунд типичным значением шага сетки по z в окрестности границы z = 0 было hz = 0.07 мкм, при этом область мелкого шага имела толщину порядка 30 мкм. В области, занимаемой платиной, толщиной порядка 10 мкм около границы z = 0 имеем hz = 0.25 мкм. В области ударной волны шаги по r и z были порядка нескольких микрометров, а в остальных областях, удаленных от поверхности мишени, — более 10 мкм. С течением времени шаги сетки постепенно увеличивались, и к концу расчета сетка становилась равномерной с шагом 20 мкм.

Точность расчета контролировалась использованием последовательности сеток и контролем закона сохранения энергии (на временах, когда волна не достигла границ

расчетной области). Заметим, что используемая схема дивергентна относительно вычисления плотности и по отношению к членам ^у (Уи - Х^Т) в уравнении (14), однако закон сохранения энергии точно не выполнялся. В приводимых ниже результатах расчетов погрешность выполнения закона сохранения энергии не превышала 0.5 %. Анализ результатов позволяет утверждать, что погрешность приводимых ниже расчетов не превышала 1-2 %.

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

6. Результаты расчетов

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

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

На рис. 1 представлены распределения плотности, давления и температуры в момент времени 47 нс (на этом и других рисунках первоначальная плотность газа соответствует 1, температура представлена в градусах, давление — в атмосферах, скорость — в километрах на секунду, расстояние — в микрометрах). На этих временах размер области высокой температуры и образующейся области низкой плотности еще мало отличается от начального размера и составляет примерно 100 мкм.

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

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

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

К моменту времени £ = 650 нс в окрестности начала координат сформирована область высокой температуры и низкой плотности (каверна) радиусом порядка 400 мкм,

Рис. 1. Распределения плотности, давления и температуры аргона в момент времени t = 47 нс при Prc = 0

0 -0.1

* -5С)(!\

-2504 ,^400 0 800 ''

Рис. 2. Распределение ¿-компоненты скорости аргона в момент времени Ь = 650 нс при Ргс = 0

размер которой в дальнейшем практически не увеличивается со временем. Плотность газа в каверне достигает своего минимального значения (в 30 раз меньше, чем первоначальная плотность аргона) в момент времени t ^ 1 мкс. Ударная волна удаляется от каверны и постепенно ослабевает (рис. 3).

Рассмотрим эволюцию распределения температуры в платине. Коэффициент поглощения платины таков, что при t = 0 вся избыточная энергия в платине сосредоточена в тонком слое у ее поверхности и температура поверхности соответствует газообразной фазе. Большие градиенты температуры вдоль z приводят к быстрому отводу тепла в глубь металла. Максимум температуры уменьшается до значений, соответствующих твердой платине, за время, равное примерно 2 нс. Температура аргона на этот момент составляет десятки тысяч градусов. Заметим, что поведение максимальной температуры поверхности платины практически одинаково для Prc = 0 и Prc = 0.1. В настоящей работе мы не рассматриваем процессы, происходящие в газообразной фазе платины. Толщина испаренного слоя платины, нагретой до температуры кипения, не превышает 0.1 мкм. Такое количество испаренного материала не окажет заметного влияния на рассматриваемую здесь динамику фонового газа, вызванную лазерным пробоем. В свою очередь динамика фонового газа может существенным образом повлиять на движение продуктов испарения, вызвав их частичное перенапыление на мишень, что составляет предмет будущих исследований.

Подчеркнем, что на описываемом этапе развития течения через кольцевую поверхность соприкосновения ударной волны с платиновой мишенью передается значительная (а на некоторых временах — наибольшая) часть тепла от аргона к платине. Для иллюстрации на рис. 4 показано распределение величины qr вдоль граница аргон—платина в различные моменты времени. Здесь q = —xdT/ dz — поток тепла в безразмерном виде,

Рис. 3. Распределения плотности, давления и температуры аргона в момент времени £ = 4.7 мкс при Ргс = 0

р, атм 1.2

1

0 г

вычисленный с учетом разрыва коэффициента теплопроводности при переходе через границу аргон—платина.

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

На временах, больших нескольких десятков микросекунд, наблюдается медленное расширение каверны и области повышенной температуры до размеров расчетной обла-

J t - 250 не A t = 650 не /1 i = 1770 не

0 400 800 1200 г

Рис. 4. Распределение потока тепла на границе аргон—платина: qr = -xrdT/dz по радиусу в различные моменты времени при Prc = 0

Рис. 5. Зависимость тепловой Цд и кинетической энергий аргона от времени при Ргс = 0 и 0.1

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

На временах > 0.1 с, когда температура аргона приближается к ее начальному значению (20 °С), нагрев платины сменяется ее охлаждением и накопленная тепловая энергия платины начинает постепенно уменьшаться.

Рассмотрим энергетические характеристики процесса. На рис. 5 представлены зависимости от времени кинетической энергии Ц и внутренней энергии газа Ц (за вычетом начальной энергии газа, соответствующей температуре Т0), нормированных на полную энергию лазерного излучения. В первые моменты времени тепловая энергия расходует-

0.45

0.35-■-■-■-■-■-——■- ,

0 2 -107 4 -107 6 ■ 107 8 ■ 107

не

Рис. 6. Зависимость тепловой энергии платины от времени при Ргс = 0 и 0.1

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

На рис. 6 показано поведение тепловой энергии платины (за вычетом начальной энергии, соответствующей температуре T0), нормированной на полную энергию лазерного излучения при разных значениях энерговыделения в газе (параметра Prc). Видно, что зависимость от Prc существенна. При Prc = 0.1 выделившаяся в платине энергия составляет 60 % энергии излучения, что соответствует экспериментам [5-7]. Это указывает на то, что эффект интерференции падающего и отраженного излучения играет заметную роль в процессе пробоя фонового газа.

Заключение

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

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

Список литературы

[1] Chichkoy B.N., Momma C., Nolte S. et al. Femtosecond, picosecond and nanosecond laser ablation of solids // Appl. Phys. A. 1996. Vol. 63. P. 109-115.

[2] Margetic V., Pakuley A., Stockhaus A. et al. A comparison of nanosecond and femtosecond laser-induced plasma spectroscopy of brass samples // Spectrochim. Acta B. 2000. Vol. 55. P. 1771-1785.

[3] Le Harzic R., Huot N., Audouard E. et al. Comparison of heat-affected zones due to nanosecond and femtosecond laser pulses using transmission electronic microscopy // Appl. Phys. Lett. 2002. Vol. 80. P. 3886-3888.

[4] Feng Q., Picard Y.N., Liu H., Yalisoye S.M., Mourou G., Pollock T.M. Femtosecond laser micromachining of a single-crystal superalloy // Scripta Material. 2005. Vol. 53. P. 511-516.

[5] Vorobyev A.Y., Guo C. Direct observation of the enhanced residual thermal energy coupling to solids in femtosecond laser ablation // Appl. Phys. Lett. 2005. Vol. 86. P. 0119161-011916-3.

[6] Vorobyev A.Y., Guo C. Enhanced energy coupling in femtosecond laser-metal interactions at high intensities // Opt. Express. 2006. Vol. 14. P. 13113-13119.

[7] Vorobyev A.Y., Kuzmichev V.M., Kokody N.G. et al. Residual thermal effects in Al following single ns- and fs-laser pulse ablation // Appl. Phys. A. 2006. Vol. 82. P. 357-362.

[8] Барчуков А.И., Бункин Ф.В., Конов В.И., Прохоров А.М. Низкопороговый пробой воздуха вблизи мишени излучением С02-лазера и связанный с ним высокий импульс отдачи // Письма в ЖЭТФ. 1973. Т. 17, № 8. С. 413-415.

[9] Pirri A.N., Root R.G., Wu P.K.S. Plasma energy transfer to metal surfaces irradiated by pulsed lasers // AIAA J. 1978. Vol. 16. P. 1296-1304.

[10] McKay J.A., Bleach R.D., Nagel D.J. et al. Pulsed-C02-laser interaction with aluminium in air: Thermal response and plasma characteristics //J. Appl. Phys. 1979. Vol. 50. P. 3231-3240.

[11] Кавашин А.В., Никитин П.И., Марине В., Сентис М.Л. Электрические поля лазерной плазмы при оптическом пробое воздуха вблизи различных мишеней // Квантовая электроника. 1998. Т. 25, № 1. С. 26-30.

[12] Агеев В.П., Барчуков А.И., Бункин Ф.В. и др. Нагрев металлов излучением импульсного СО2-лазера // Квантовая электроника. 1979. Т. 6, № 1. С. 78-85.

[13] Mao S.S., Mao X., Greif R., Russo R.E. Dynamics of an air breakdown plasma on a solid surface during picosecond laser ablation // Appl. Phys. Lett. 2000. Vol. 76. P. 31-33.

[14] Климентов С.М., КононЕнко Т.В., Пивоваров П.О. и др. Роль плазмы в абляции материалов ультракороткими лазерными импульсами // Квантовая электроника. 2001. Т. 31, № 5. С. 378-382.

[15] Ireland C.L.M., Morgan C.G. Gas breakdown by a short laser pulse // J. Phys. D: Appl. Phys. 1973. Vol. 6. P. 720-729.

[16] Bulgakova N.M., Zhukov V.P., Vorobyev A.Y., Guo Ch. Modeling of residual thermal effect in femtosecond laser ablation of metals: role of gas environment // Appl. Phys. A. 2008. in press, DOI 10.1007/s00339-008-4568-1.

[17] Федорук М.П., Березин Ю.А. Моделирование нестационарных плазменных процессов. Новосибирск: Наука, 1993.

[18] Калиткин Н.Н., Альшин А.Б., Альшина Е.А., Рогов Б.В. Вычисления на квазиравномерных сетках. М.: Физматлит, 2005.

[19] Жуков В.П. Конечно-разностная схема для решения двухжидкостных МГД-уравнений в цилиндрической системе координат // Журн. вычисл. математики и мат. физики. 2005. Т. 45, № 1. С. 156-169.

Поступила в редакцию 15 февраля 2008 г., в переработанном виде — 6 марта 2008 г.