Сравнительный анализ оптоакустических импульсов в алюминии и кремнии
УДК 538.975: 548.25
Сравнительный анализ оптоакустических импульсов в алюминии и кремнии
© 2011 г. О. Н. Королева*, канд. физ.-мат. наук, А. В. Мажукин** ** Московский гуманитарный университет, Москва ** Институт прикладной математики им. М.В. Келдыша РАН, Москва ** Е-mail: koroleva.on@mail.ru, specimen@modhef.ru
Методами математического моделирования проведен сравнительный анализ о птоакустического отклика на короткоимпульсное лазерное воздействие умеренной интенсивности на алюминий и кремний. Рассмотрены режимы воздействия, близкие к экспериментальным данным, когда использовались лазерные импульсы одинаковой длительности 3×10–9 с, интенсивность которых менялась в пределах 107–1010 Вт/см2.
Ключевые слова: математическое моделирование, оптоакустика, лазерное воздействие, динамическая адаптация.
Коды OCIS: 350.3390, 320.4240, 350.5340
Поступила в редакцию 18.02.2011
Введение
Бурное развитие лазерных технологий в различных областях обработки и создания новых материалов в последние десятилетия привлекает внимание к исследованию процессов, протекающих в зоне лазерного воздействия и к основному средству такого исследования – оптоакустической диагностике [1–4]. Опто акустическая диагностика является точным и чувствительным способом всестороннего описания динамики протекающих процессов. Наибольший интерес для лазерных технологий представляет импульсное лазерное воздействие умеренной интенсивности (G ≤ 1010 Вт/см2). Действие лазерного излучения на поглощающие конденсированные среды сопровождается генерацией импульсов давления, несущих информацию о характере процессов в зоне облучения. В основе оптоакустики лежит измерение акустических полей, возбуждаемых лазерным воздействием в конденсированных средах. Измерение акустических сигналов в рамках одномерной по пространству модели сводится к регистрации временной развертки импульса давления P(t). Благодаря относительной простоте экспериментальной реализации акустические методы перекрывают широкий диапазон ин тенсивности 106–1010 Вт/см2 и длительности
воздействия 10–9–10–3 с. В ряде прикладных задач акустические методы позволяют получать информацию о процессах в конденси- рованных средах, недоступную для других способов измерений. Оказалось, что изучение формы акустических сигналов является эффективным способом исследования динамики быстрых фазовых переходов в процессе воз действия нано- и микросекундных лазерных импульсов на конденсированные среды [5].
При интерпретации экспериментально полученных сигналов трудно установить сочетанием каких процессов они были образованы, если отсутствует качественная информация о поведении и относительном вкладе в сигнал тех или иных процессов. Инструментом полу чения такой информации являются современные методы математического моделирования.
Работа посвящена изучению особенностей оптоакустических сигналов, возникающих в связи с плавлением и испарением вещества в зоне лазерного облучения алюминия и кремния.
Постановка задачи
На свободную поверхность плоской сильнопоглощающей пластины воздействует лазерный импульс наносекундной длительности с интенсивностью 107–109 Вт/см2. Энергия
“Оптический журнал”, 78, 8, 2011
79
vsl vsur
x = x0
Твердая фаза
Гsl
Жидкая фаза
Гlv
Лазерное излучение
Рис. 1. Пространственное расположение фаз и межфазных границ при лазерном воздействии
на мишень: Гlv – граница испарения; Гsl – граница раздела твердое тело–жидкость; vsl, vlv – скорости движения границ испарения и плав-
ления.
импульса поглощается в тонкой приповерхностной зоне мишени, что вызывает ее нагрев, плавление и испарение. Под действием теплового расширения и фазовых превращений вблизи поверхности формируется сигнал давления (оптоакустический сигнал), который распространяется к противоположной закрепленной сторонепластины.Радиуспятнафокусировкиrлазерного излучения предполагается много большим глубины проникновения лазерного излучения r >> (χτL)1/2, где χ – коэффициент температуропроводности (для кремния χ = 1,01 см2/с; для алюминия χ = 0,99 см2/с), τL – полуширина лазерного импульса. Распределение интенсивности в плоскости пятна предполагается равномерным, что позволяет использовать
для математического описания одномерное по пространству приближение (рис. 1).
Математическое моделирование процессов формирования и распространения оптоакустического сигнала осуществляли в рамках совмещенного варианта гидродинамической задачи Стефана [6], включающего классический и однофазный вариант. Основу совмещенного варианта составляет система уравнений гидродинамики, дополненная уравнением переноса лазерного излучения в материале. Cистема гидродинамики включает уравнения неразрыв ности, движения и полной энергии. Учитываются конвективный, кондуктивный и излучательный механизмы переноса. В дивергентной форме система уравнений имеет вид
ëêêé
¶ρ ¶t
+
¶ ¶x
(ρu)=
0ùúúûk,
( )êëéê
¶ ¶t
(ρu)+
¶ ¶x
ρu2
=
-
¶p ¶x
úûùú
k
,
(1) (2)
ëêêêé
¶ ¶t
èççççæρèççççæH
+
u2 2
ø÷÷÷÷öø÷÷÷÷÷ö
+
¶ ¶x
èççççæρuèççççæH
+
u2 2
ø÷÷÷ö÷ø÷÷÷÷÷ö
=
=
-
¶ ¶x
(
pu)-
¶W ¶x
-
¶G ¶x
ùûúú
k
,
(3)
ëêêé
¶G ¶x
+
α
G
=
0ûúúùk
x Î[x0,Ãsl ][Ãsl, Ãlv ]
k = ïíïïìïîsl,, xxÎÎ[[Ãxs0l,,ÃÃslvl ]].
В качестве уравнения состояния использовали соотношения
Hl = cplTl
pl
=
uc2ρ0l
ççæçèççæççè
ρ ρ0l
-1÷÷öø÷÷+ βl
(T
-Tm
)÷÷ö÷÷ø
x Î[Ãlv, Ãsl ].
Hs = cpsTs
ps
=
uc2ρ0s
èççççæèçççæ
ρ ρ0s
-1ø÷÷÷ö÷
+
βs
(T
-T0
)÷ø÷÷÷ö
x Î[Ãsl, xs ].
(4)
(5) (6)
В рассматриваемой системе уравнений
использованы следующие обозначения и единицы: ρ [г/см3] – плотность, u [см/с],
р [бар] – гидродинамическая скорость и давление, uс[см/с] – скорость звука, H [Дж/см3] – энтальпия, W [Вт/см2] – тепловой поток,
α [1/см] – коэффициент поглощения лазерного излучения, G [Вт/см2] – плотность мощно-
сти лазерного излучения, cp [Дж/г K] – удельная теплоемкость при постоянном давлении,
β [1/K] – температурный коэффициент линей-
ного расширения, T0 [K] – комнатная темпе ратура, Tm [K] – равновесная температура плавления, ρ0, s, ρ0, l – начальные плотности твердой и жидкой фазы, задаваемые соответ-
ственно при температурах T0 и Tm. Индексы sur, s, l, υ обозначают принадлежность ве-
личин соответственно к поверхности, твердой,
жидкой и парообразной средам, а индексы sl
и lυ – принадлежность к межфазным грани-
цам раздела твердое тело–жидкость x = Γsl и границе испарения Гlυ. Система уравнений
80 “Оптический журнал”, 78, 8, 2011
(1)–(6) дополнялась граничными и начальными
условиями.
Левая граница х = х0, полагалась закрепленной и теплоизолированной (рис. 1).
x = x0:
λ
¶Ts ¶x
= 0,
us
= 0.
(7)
На межфазной границе x = Γsl формулируется модель поверхностного плавления, состоя-
щая из трех законов сохранения, дополненных
условием непрерывности температуры фазово-
го перехода.
x = Ãsl: ρsusl = ρl (us -ul + usl ), ρsu2sl + ps = ρl (us -ul + usl )2 + pl,
(8) (9)
λl
¶Tl ¶x
- λs
¶Ts ¶x
=
Lnme ρs usl ,
(10)
Ts = Tl = Tsl = Tm,
(11)
где
( )( )Lnme = Lm + cpl - cps
Tsl -Tm, 0
+
ρs ρs
+ -
ρl ρl
(us
- ul )2
2
– неравновесная теплота плавления, Lm [Дж/г] – равновесная теплота плавления, λ [Вт/см K] –
коэффициент теплопроводности. На правой облучаемой границе x = Γlυ(t) фор-
мулируется модель поверхностного испарения,
записанная в приближении кнудсеновского
слоя и состоящая из трех законов сохранения
и двух дополнительных условий, характери
зующих степень неравновесности фазового пе-
рехода:
x = Ãlu : ρlulu = ρu (ulu + ul -uu ),
(12)
pl + ρlu2lu = pu + ρu (ulu + ul -uu )2 ,
(13)
WlT = λl
¶Tl ¶x
= ρluluLnue
+ σTs4ur ,
(14)
Tu = TlαΤ (M), ρu = ρsatαρ (M),
(15)
Glu
=
A(Tsur
)G0
expèççççæç-4çèçæç
t tL
ø÷÷÷÷ö2
÷ø÷÷÷÷÷ö,
(16)
где
Lnue
=
Lu
(Tl
)+
cpu
(Tl
-Tu)+
ρl ρl
+ -
ρu ρu
(ul
- uu)2
2
– неравновесная теплота испарения, Glυ, G0 – поглощенная поверхностью энергия лазер-
ного излучения и ее максимальное значение,
Asur [%] – поглощательная способность поверхности, τL [нс] – полуширина лазерного импульса, Lυ [Дж/г] – теплота кипения, M = = ulυ/uc – число Маха, uc = (γRTυ)1/2 – скорость звука на внешней стороне кнудсеновского слоя,
γ = cp/cυ, R [Дж/г K] – универсальная газовая постоянная, αT(M), αρ(M) – коэффициенты Крута [7, 8], при М = 1 αT(M) = 0,633, αρ(M) = 0,326.
Давление насыщенного пара вычисляли, ис-
пользуя уравнения Клапейрона–Клаузиуса:
psat
=
pb
expæççççè
Luµ RTb
æçççè1-
Tb Tl
ö÷÷÷ø÷÷ö÷÷÷ø,
ρsat = psat /RTl,
(17)
где pb, Tb – соответственно давление и температура кипения при нормальных условиях. В на-
чальный момент времени t0 температуру и плотность полагали постоянными, а скорость,
равной нулю:
Ts = T0, ρs = ρ0,s, us = 0.
(18)
Оптические и теплофизические характеристики алюминия и кремния
На форму оптоакустических сигналов большое влияние оказывают теплофизические и оптические характеристики материала мишени [9, 10]. Используемые в расчетах характеристики материалов, представленные на рис. 2a–б и 3, взяты из справочных данных [11, 12]. Вертикальными линиями на рисунках отмечены равновесные температуры плавления и испарения каждого из материалов. Все теплофизические и оптические характеристики обоих материалов претерпевают разрыв при переходе через значение равновесной температуры плавления (на рис. 2, 3 она отмечена вертикальными линиями). Учитывая, что один из рассматриваемых материалов является металлом, а другой полупроводником, их теплофизические и оптические характеристики сильно различаются. При переходе через фронт плавления теплофизические характеристики алюминия скачком уменьшаются, а кремния – возрастают. В таблице приведены значения теплофи зических параметров, характеризующие свойства рассматриваемых материалов.
Алгоритм решения. Основная сложность решения задачи Стефана заключается в наличии двух подвижных границ, положение которых неизвестно и должно определяться в ходе расчетов, при этом размеры твердой и жидкой подобластей в процессе расчетов могут меняться на несколько порядков. Анализ вклада
“Оптический журнал”, 78, 8, 2011
81
, Вт/см3 , г/см K
3
2
1
Tm
1000
, Вт/см3 , г/см K
(а)
1
Tu
3000
5000
(б)
2 3
7000
Ср, Дж/гK
3 2 1
Ср, Дж/гK
31
3
22
3 2
11
1000
Tm
2000
Tu 3000 4000 Т, K
Рис. 2. Теплофизические характеристики алюминия (а) и кремния (б). 1 – плотность; 2 – те-
плоемкость; 3 – теплопроводность; Тm – температура плавления; Тν – температура испа рения.
фазового перехода в оптоакустический сигнал требует точного определения положения межфазной границы и ее скорости. Наиболее эффективным методом для решения поставленной задачи является метод динамической адаптации [13–16], позволяющий производить расчет с явным выделением межфазных границ.
Метод динамической адаптации. Численное решение задачи (1)–(18) осуществляется методом динамической адаптации, в основу которого положена идея перехода к произвольной нестационарной системе координат, осуществляемого с помощью искомого решения. В про-
100
A, %
80
60
(а)
40
20
0
Tm
1000
Tu
3000
5000
7000
100
A, %
80
60
(б)
40
20
Tm
0
1000
2000
3000 Т, K
Рис. 3. Поглощательная способность алюминия (а) и кремния (б). Тm – температура плавления. Тv – температура испарения.
извольной нестационарной системе координат проблема описывается расширенной дифференциальной системой уравнений, часть из которых описывает физическое явление, а другая часть – динамику узлов расчетной сетки. Проблемы, связанные с подвижными границами, снимаются путем перехода к произвольной нестационарной системе координат, в которой узлы сетки и границы оказываются неподвижными. Преобразование координат осуществляется автоматически с помощью искомого решения, что позволяет производить размещение узлов сетки в зависимости от особенностей решения.
Переход из физического пространства с переменными (x, t) в расчетное с произвольной нестационарной системой координат (q, τ) осу-
Теплофизические параметры алюминия и кремния
Минерал
Атомный вес, г/моль
Al 26,98 Si 28
Температура
Tm K
933 1683
Tb K
2793 3514
Коэффициент теплового расширения, 1/K
Твердая фаза βs
2,33×10–5 4,65×10–3
Жидкая фаза βl
5,00×10–5 1,40×10–4
Теплота перехода
Плавление Lm (Дж/г)
400,30 1797,0
Кипение Lb (Дж/г)
10860 13720
82 “Оптический журнал”, 78, 8, 2011
ществляется с помощью замены переменных общего вида
x = f(q, t), t = t,
(19)
имеющей обратное невырожденное преобразо-
вание q = j(x, t), t = t. Частные производные в
системе координат (q, τ) имеют следующий вид:
¶ ¶x
=
1 ψ
¶ ¶q
,
¶ ¶t
=
¶ ¶t
+
Q ψ
¶ ¶q
,
(20)
где ∂x/∂τ = –Q – скорость движения новой си-
стемы координат относительно исходной, под-
лежащая в дальнейшем определению.
В нестационарной системе координат две
ф азовые подобласти [Гlν, Гsl] ∪ [Гsl, xs] с по движными границами Гlν, Гsl отображаются на две расчетные подобласти [0, qsl] и [qsl, qs] с неподвижными границами. До появления
жидкой фазы физическая область [Гlν, xs] отображается на расчетную подобласть [0, qs]. Физическая координата х в расчетном простран-
стве становится новой неизвестной функцией.
В новой системе координат динамика расчет-
ной сетки описывается дополнительным диф-
ференциальным уравнением:
êëéê
¶ψ ¶t
=
-
¶Q ¶q
úúùû k
k = s, l,
(21)
где ψ = ∂x/∂q – метрический коэффициент, Q – функция, определяющая конкретный вид
преобразования в соответствии с особенностями
задачи.
Используя преобразование переменных (19),
запишем дифференциальную задачу (1)–(17)
в произвольной нестационарной системе координат (q, τ):
êëêé
¶(ψρ)
¶t
+
¶ ¶q
(ρ(u
+
Q))=
0ùûúú
,
k
êëêé
¶ ¶t
(ψρu)+
¶ ¶q
(ρu(u
+
Q))=
-
¶p ¶q
ùúúû
,
k
éêêêë
¶ ¶t
æçççèçψρæççççèH
+
u2 2
ö÷ø÷÷÷ö÷ø÷÷÷÷
+
¶ ¶q
æçççèçρ(u
+
Q)æççççè
H
+
u2 2
÷ö÷÷÷ø÷÷÷÷÷öø
=
=
-
¶ ¶q
(
pu)+
¶ ¶q
ççèçæç
λ(T)
ψ
¶T ¶q
÷ø÷ö÷÷ûúúúù
,
k
(22)
êéëê
¶G ¶q
+
αψG
=
0úùúû
,
k
q Î[0, qsl ][qsl, qs ], t Î[t0, tend ], k = s, l.
Граничные условия. q = q0: u = W = 0, Q = 0;
(23)
q = Γsl: ρs (Qsl + us )= ρl (Qsl + ul ), ps + ρs (Qsl + us )2 = pl + ρl (Qsl + ul )2,
(24) (25)
Wl - Ws = ρsuslLnme, usl = Qsl + us,
q = Γlu: ρl (Qlu + ul )= ρu (Qlu + uu ), pl + ρl (Qlu + ul )2 = pu + ρu (Qlu + uu )2,
(26) (27) (28)
WlT
=
WlT
-
λl ψl
¶Tl ¶q
= ρluluLnue + σTs4ur ,
ulu = -(Qlu + ul ),
Glu
=
A(Tsur
)G0expçèçççæç-4èçççæ
t tL
÷÷ø÷÷ö2
÷÷ø÷÷ö÷÷.
(29)
t = t0: Ts = T0, us = 0, ρs = ρ0, s, ψ = 1. (30)
Таким образом, при переходе к произвольной нестационарной системе координат исходная дифференциальная модель трансформируется в расширенную дифференциальную систему (21)–(30), в которой появляется дополнительное уравнение типа (21), являющееся уравнением обратного преобразования. Его тип, свойства и вид краевых условий зависят от конкретного вида функции Q. Для построения равномерных (квазиравномерных) на каждый момент времени сеток в областях с подвижными границами функция Q задается в виде [13, 14]:
êéêëQ
=
-D
¶ψ ¶q
ûùúú
,
k
k = s, l,
где коэффициент диффузии D выражается че-
рез геометрические и скоростные параметры
задачи: D = L2(t)(|υsl| + |υlυ|)/Δx. Разностная аппроксимация. Численное ре-
шение дифференциальной модели (21)–(30)
осуществлялось при помощи конечно-разност-
ного метода, согласно которому уравнения ги-
дродинамики аппроксимировали семейством
консервативных разностных схем, получен-
ных интегро-интерполяционным методом [17].
Для построения семейства разностных схем в расчетном пространстве вводится сетка wqt с неравномерными шагами hk, i, τj по пространственной q и временной τ переменным:
ω = {ωl ∪ ωs}×{ωτ},
“Оптический журнал”, 78, 8, 2011
83
где
wl
=
îïïìïíïiq=l,
i, ql, i+1/2; 0,.., Nl -1,
ql, i+1 = ql, ql, 0 = 0,
i+ qNl
hl, =
i+1, qsl ,
ql, i+1/2 hl, 0 = 0,
= ql, i hl, Nl
+
+1
0,5hl, =0
i+1
ïüïþïýï,
ws
=
ìïïïîïíiq=s,
i, 0,
qs, ..,
i+1/2; Ns -1,
qs, i+1 = qs, i + hs, i+1, qs, 0 = 0, qNs = qsl,
qs, i+1/2 hs, 0 = 0,
= qs, i hs, Ns
+
+1
0,5hs, =0
i+1
ýþüïïïï,
{ }wt = tj, tj+1 = tj + Δtj+1, j = 0,..., J, t0 = t0 .
При разностной аппроксимации уравнений
к узлам с целочисленными индексами отно
сятся сеточные функции p, Q, x, к точкам с полуцелыми индексами (qm+1/2, τj) – сеточные функции ψ, ρ, u, T. Аппроксимация системы
уравнений (22) производилась с помощью полностью неявной разностной схемы, имеющей порядок аппроксимации Q(τ + h).
Для дифференциальной системы (21), (22) семейство разностных схем имеет вид:
ïíïîïïïïì(ψρ)ij++11/Δ2t-j+(1ψρ)ij+1/2
+
ëéρ(u
+
Q)ûùij++11 - ëéρ(u
hi+1
+
Q)ûù
j+1 i
=
0ýïïïïïþüïk,
(31)
ïíîïïïïïì(ψρu)jj++11/Δ2t-j+(1ψρu)jj+1/2
+
ëéρu(u
+
Q)ûù
j+1 i+1
-
ëéρu(u
+
Q)ûù
j+1 i
hi+1
=-
pij++11 + pij+1 hi+1
ïþïïïýïïük,
(32)
( ) ( ) ( ) ( )ïíïïïïïïïïîïïïïïïïì=éêëψ-ρ
ε + u2
2
ùúû
j+1 i+1/2
-
éëêψρ
ε
Δtj+1
(pu)ij++11 +(pu)ij+1 + 1
hi+1
hi
+ u2
2
ùûú
j i+1/2
+ éêêëρ(u + Q)
ε + u2
2
ij++11ùúúû - éêëêρ(u
hi+1
( )ççççèçæ
λij++11
ψ
j+1 i+1
Tij++31/2 -Tij++11/2 0,5 hi+1 + hi
-
λij+1 ψij+1
T0i,j+5+1(1/h2i-+Thi-ji+-111/2)ø÷÷÷÷÷ö
+
Q)
ε + u2
2
ij+1ùúúû
=ïïýïïïïïïïïïïüïþïïïk,
(33)
éêëêê
ψij++11/2 - ψij+1/2 Δtj+1
=
-
Qij++11 - Qij+1 hi+1
ùúûúú
,
k
(34)
ëêêêé
Gij++11 - Gij+1 hi+1
+
α(ψG
)j+1
i+1/2
=
0úúûúù
,
k
(35)
éêêêëQij+1
=
-Dij
ψ
j+1 i+1/2
-
ψij-+11/2
(hi + hi+1 ) 2
úùúûú,
ëêêêé
xij++11 - xij+1 hi+1
=
ψij++11/2
ûúùúú
,
k
где i = 0, …, Nl – 1 при k = l, i = 0,…,Ns – 1 при k = s.
Аппроксимация граничных условий (23)–
(30) выполнялась следующим образом:
84
q = q0:
ψsj,+01 λsj,+01
Tsj,+01 -Tsj,+1/12 hs, 0 /2
= 0,
us, 0 = 0, Qs0 = 0,
(36)
q = Γsl : (ρsVsl )j+1 = ρlj+1 (us - ul + vsl )j+1, (37)
( )ρsVs2l j+1+ psj+1=
( )= ρlj+1
(us -ul + vsl )j+1
2
+
plj+1.
(38)
λs, N-1 ψs, N-1
Tsj,+N1-1/2 -Tsj,+N1-1 hs, N-1/2
-
λl, 0 ψl, 0
Tl,j+01 -Tl,j1+/12 hl,1/2
=
= Lm (ρsvsl )j+1,
(39)
Ts = Tl = Tm.
“Оптический журнал”, 78, 8, 2011
q = Γlu:
(ρvvlv )j+1= (ρv (ulv -uv + vlv ))j+1,
( )Qljv+1 = - vljv+1 + uljv+1 ,
(40)
( )ρlvvl2v j+1+ pljv+1=
( )= ρv (ulv -uv + vlv )2 j+1+ pvj+1,
(41)
λsj+ur1
( )ψsj+ur1
Tkj,+11/2 -Tsju+r1 hk,1/2
= (ρsur LvVkv )j+1+
σT4
j+1,
(42)
( )(Glu )j+1= A Tsju+r1
G0
expçèççæçç-4çèæçç
t tL
÷ø÷÷÷ö2
÷÷ø÷÷÷÷öj+1.
(43)
Величины в полуцелых узлах для сеточных функций, отнесенных к целым точкам, вычисляются как полусумма значений в ближайших узлах. При определении величин в узлах с целочисленными индексами для функций, о тнесенных к полуцелым точкам, используется линейная интерполяция:
ym
=
ym-1/2ψm+1/2 ψm+1/2
+ ym+1/2ψm-1/2 + ψm-1/2
.
(44)
Результаты моделирования
Рассмотрим и проанализируем изменения амплитудных характеристик оптоакустических сигналов в алюминии и кремнии, вызванных короткоимпульсным лазерным воздействием.
Моделирование формы акустического си гнала проводилось для лазерного импульса гауссовой формы с длительностью на полувысоте τL = 3 нс, в диапазоне пиковых интенсивностей G0 = 107 – 1010 Вт/см2. Кривые давления получены для глубины хℓ =10 мкм. Точка расположена вне зоны теплового влияния импульса и на рассматриваемом интервале времени основной сигнал в ней не искажается отраженным. Выбранный диапазон интенсивности излучения для данной длительности τL соответствует области доминирования процесса плавления обоих материалов, а испарение еще не вносит существенного вклада. Сравнительный анализ оптоакустических сигналов должен дать представление:
–– о пороге плавления обоих материалов; –– о характерном профиле оптоакустических сигналов для металлов и полупроводников; –– о влиянии теплофизических и оптических параметров облучаемых материалов на динами-
ку фазовых переходов и амплитуды оптоакустических сигналов.
Для теплового механизма генерации звука существуют оценки [1], полученные из совместного решения линейных уравнений теплопроводности и гидродинамики, согласно которым:
P(t)=
χβρ0
¶Tsur ¶t
+
1 cpα
¶G¶t(t),
χ
=
λ cpρ0
.
(45)
На рис. 4а–в представлены временные за-
висимости температуры поверхности, скоро-
сти фронта плавления и давления для крем-
ния. Моделирование показало, что, несмотря
на более высокую, чем у алюминия, темпе
ратуру плавления TmSi (1683 K)> TmAl(933,6 K),
порог плавления кремния оказался ниже
(рис. 4a). Пороговая интенсивность лазерного воздействия составила G0 = 2,4×107 Вт/см2. Температура плавления на облучаемой поверх-
ности кремния достигается на заднем фронте
лазерного импульса (рис. 4a). Быстрому нагре-
ву приповерхностных слоев кремния способству-
ет значительное уменьшение теплопроводности
в твердой фазе с ростом температуры. Появле-
ние жидкой фазы характеризуется скоростью Vsl (рис. 4б), отрицательная ветвь которой характеризует плавление, а положительная – за-
твердевание. Оптоакустический сигнал явля-
ется динамической характеристикой, которая
весьма чувствительна к темпу ввода энергии в
систему. На начальной стадии нагрева импульс
давления (рис. 4в), положителен. Величина его
растет с ростом температуры до значения рa ≈ 3,5 бар, а затем замедляется из-за уменьшения производных ∂Tsur/∂t, ∂G(t)/∂t и плавно спадает после смены их знака. Момент обра-
зования жидкой фазы характеризуется отрицательным скачком плотности Δρ = ρs – ρl 0, отрицательная полуволна – сжатию ∂Ts/∂t 0. Это проявляется на развертке давления в виде резкого коротко-
го скачка вверх, вплоть до нулевого значения давления рa. Последующее затвердевание происходит при быстро уменьшающемся энергов-
кладе, в результате чего максимальная скорость
затвердевания (≈ 60 м/с) в 1,5 раза превышает
максимальную скорость плавления (рис. 5б), и в
7,5 раз превышает максимальную скорость за-
твердевания кремния. Столь высокая скорость
затвердевания вызывает резкое увеличение отрицательного давления, рa ≈ 15 бар (рис. 5в).
При повышении интенсивности излучения до G0 = 2×108 Вт/см2 оптоакустические сигна-
86 “Оптический журнал”, 78, 8, 2011
Тsur, K
800
400
0 – 4 – 2
0
(а) usl, м/с
(б) Р, бар
40 10
20 0
0
– 20 – 10
– 40 – 20
2 4 6 t, нс
– 4 – 3 – 2 –1 0 1 2 3 t, нс
– 4 – 2 0
2
(в) 4 6 t, нс
Рис. 5. Зависимость от времени температуры поверхности Тsur, скорости плавления usl и давления р для алюминия. Пунктиром обозначена зависимость поглощенной энергии лазерного излучения от времени. Полуширина лазерного импульса 3 нс. G0 = 3,5×107 Вт/см2.
Тsur, K
2500 2000
1500
(а)
usl, м/с
300
200
100
0
1000
– 100
500 – 200
0 – 2
0
2
4
– 300
6 t, нс
– 4
04
(б) ulv, м/с
8 10–3 4 10–3 0
– 4 10–3 – 8 10–3 – 1,2 10–2
8 t, нс
Р, бар
200 100
0 – 100
– 2
0
(в) 2 4 t, нс
Рис. 6. Зависимость от времени температуры поверхности Тsur, скорости плавления usl, испарения ulν и давления р для алюминия. Пунктиром обозначена зависимость поглощенной энергии от времени. Полуширина лазерного импульса 3 нс. G0 = 2×108 Вт/см2.
Тsur, K
2400
1600
800
0 – 4 – 2
0
(а)
24
usl, м/с
2 10–3
(б)
0
– 2 10–3
– 4 10–3
– 6 10–3
– 8 10–3
– 1 10–2
6 t, нс
– 4 – 2 0 2 4
ulv, м/с
4 10–5
Р,
бар
10
0 – 4 10–5
0
– 8 10–5
– 10
– 1,2 10–4
– 20
– 1,6 10–4 – 2 10–4
6 t, нс
– 30 – 40
– 4
– 2
0
2
(в) 4 t, нс
Рис. 7. Зависимость от времени температуры поверхности Тsur, скоростей плавления usl, испарения ulν и давления р для кремния. Пунктиром обозначена зависимость поглощенной энергии от времени. Полуширина лазерного импульса 3 нс. G0 = 2×108 Вт/см2.
лы алюминия и кремния качественно меняются (рис. 6a–в, 7а–в). Момент плавления перемещается на передний фронт лазерного импульса, где возрастает роль поглощательной способности материала. В алюминии поглощательная способность при плавлении скачком увеличивается, усиливая скорость плавления до 300 м/с (рис. 6б) и, как следствие, происходит скачок давления до 270 бар (рис. 6в). Поверхность при этом разогревается до температуры TsAulr ³ TbAl и становится заметным вклад испарения, скорость которого достигает vlν ≈ 8×10–3 м/с
(рис. 6б). Амплитуда оптоакустического импульса после всплеска давления, вызванного скачком поглощательной способности и положительным скачком плотности, под совместным воздействием плавления, испарения и уменьшения интенсивности излучения быстро спадает, переходя в отрицательную полуволну.
Оптоакустический сигнал в кремнии имеет еще более сложную временную структуру (рис. 7а–в). Для кремния поглощательная способность поверхности скачком уменьшается
“Оптический журнал”, 78, 8, 2011
87
(рис. 7а), что в совокупности с отрицательным скачком плотности приводит к увеличению отрицательного скачка давления до 35 бар (рис. 7в). Продолжающийся нагрев позволяет уменьшить отрицательное давление до нуля, которое затем снова увеличивается до 10 бар, из-за совместного воздействия испарения, скорость которого (vlv ≈ 4,5×10–3 м/с, рис. 7б) становится заметной, и уменьшения интенсивности лазерного импульса. Процесс затвердевания приводит к формированию положительной полуволны с максимальным значением pa ≈ 10 бар. Из-за более низкой п оглощательной способности темп нагрева кремния ниже, чем у алюминия, чем и объясняются более низкие максимальные значения температуры поверхности, скорости плавления и скачков давления.
Заключение
Выполненный с помощью математического моделирования сравнительный анализ опто акустических сигналов, вызванных лазерным воздействием на алюминий и кремний, показал, что:
1. Плавление существенно меняет вид оптоакустического сигнала в кремнии и алюминии.
2. Скачок давления, связанный с плавлением, положителен у алюминия (Δp > 0) и отрицателен у кремния (Δp
Сравнительный анализ оптоакустических импульсов в алюминии и кремнии
© 2011 г. О. Н. Королева*, канд. физ.-мат. наук, А. В. Мажукин** ** Московский гуманитарный университет, Москва ** Институт прикладной математики им. М.В. Келдыша РАН, Москва ** Е-mail: koroleva.on@mail.ru, specimen@modhef.ru
Методами математического моделирования проведен сравнительный анализ о птоакустического отклика на короткоимпульсное лазерное воздействие умеренной интенсивности на алюминий и кремний. Рассмотрены режимы воздействия, близкие к экспериментальным данным, когда использовались лазерные импульсы одинаковой длительности 3×10–9 с, интенсивность которых менялась в пределах 107–1010 Вт/см2.
Ключевые слова: математическое моделирование, оптоакустика, лазерное воздействие, динамическая адаптация.
Коды OCIS: 350.3390, 320.4240, 350.5340
Поступила в редакцию 18.02.2011
Введение
Бурное развитие лазерных технологий в различных областях обработки и создания новых материалов в последние десятилетия привлекает внимание к исследованию процессов, протекающих в зоне лазерного воздействия и к основному средству такого исследования – оптоакустической диагностике [1–4]. Опто акустическая диагностика является точным и чувствительным способом всестороннего описания динамики протекающих процессов. Наибольший интерес для лазерных технологий представляет импульсное лазерное воздействие умеренной интенсивности (G ≤ 1010 Вт/см2). Действие лазерного излучения на поглощающие конденсированные среды сопровождается генерацией импульсов давления, несущих информацию о характере процессов в зоне облучения. В основе оптоакустики лежит измерение акустических полей, возбуждаемых лазерным воздействием в конденсированных средах. Измерение акустических сигналов в рамках одномерной по пространству модели сводится к регистрации временной развертки импульса давления P(t). Благодаря относительной простоте экспериментальной реализации акустические методы перекрывают широкий диапазон ин тенсивности 106–1010 Вт/см2 и длительности
воздействия 10–9–10–3 с. В ряде прикладных задач акустические методы позволяют получать информацию о процессах в конденси- рованных средах, недоступную для других способов измерений. Оказалось, что изучение формы акустических сигналов является эффективным способом исследования динамики быстрых фазовых переходов в процессе воз действия нано- и микросекундных лазерных импульсов на конденсированные среды [5].
При интерпретации экспериментально полученных сигналов трудно установить сочетанием каких процессов они были образованы, если отсутствует качественная информация о поведении и относительном вкладе в сигнал тех или иных процессов. Инструментом полу чения такой информации являются современные методы математического моделирования.
Работа посвящена изучению особенностей оптоакустических сигналов, возникающих в связи с плавлением и испарением вещества в зоне лазерного облучения алюминия и кремния.
Постановка задачи
На свободную поверхность плоской сильнопоглощающей пластины воздействует лазерный импульс наносекундной длительности с интенсивностью 107–109 Вт/см2. Энергия
“Оптический журнал”, 78, 8, 2011
79
vsl vsur
x = x0
Твердая фаза
Гsl
Жидкая фаза
Гlv
Лазерное излучение
Рис. 1. Пространственное расположение фаз и межфазных границ при лазерном воздействии
на мишень: Гlv – граница испарения; Гsl – граница раздела твердое тело–жидкость; vsl, vlv – скорости движения границ испарения и плав-
ления.
импульса поглощается в тонкой приповерхностной зоне мишени, что вызывает ее нагрев, плавление и испарение. Под действием теплового расширения и фазовых превращений вблизи поверхности формируется сигнал давления (оптоакустический сигнал), который распространяется к противоположной закрепленной сторонепластины.Радиуспятнафокусировкиrлазерного излучения предполагается много большим глубины проникновения лазерного излучения r >> (χτL)1/2, где χ – коэффициент температуропроводности (для кремния χ = 1,01 см2/с; для алюминия χ = 0,99 см2/с), τL – полуширина лазерного импульса. Распределение интенсивности в плоскости пятна предполагается равномерным, что позволяет использовать
для математического описания одномерное по пространству приближение (рис. 1).
Математическое моделирование процессов формирования и распространения оптоакустического сигнала осуществляли в рамках совмещенного варианта гидродинамической задачи Стефана [6], включающего классический и однофазный вариант. Основу совмещенного варианта составляет система уравнений гидродинамики, дополненная уравнением переноса лазерного излучения в материале. Cистема гидродинамики включает уравнения неразрыв ности, движения и полной энергии. Учитываются конвективный, кондуктивный и излучательный механизмы переноса. В дивергентной форме система уравнений имеет вид
ëêêé
¶ρ ¶t
+
¶ ¶x
(ρu)=
0ùúúûk,
( )êëéê
¶ ¶t
(ρu)+
¶ ¶x
ρu2
=
-
¶p ¶x
úûùú
k
,
(1) (2)
ëêêêé
¶ ¶t
èççççæρèççççæH
+
u2 2
ø÷÷÷÷öø÷÷÷÷÷ö
+
¶ ¶x
èççççæρuèççççæH
+
u2 2
ø÷÷÷ö÷ø÷÷÷÷÷ö
=
=
-
¶ ¶x
(
pu)-
¶W ¶x
-
¶G ¶x
ùûúú
k
,
(3)
ëêêé
¶G ¶x
+
α
G
=
0ûúúùk
x Î[x0,Ãsl ][Ãsl, Ãlv ]
k = ïíïïìïîsl,, xxÎÎ[[Ãxs0l,,ÃÃslvl ]].
В качестве уравнения состояния использовали соотношения
Hl = cplTl
pl
=
uc2ρ0l
ççæçèççæççè
ρ ρ0l
-1÷÷öø÷÷+ βl
(T
-Tm
)÷÷ö÷÷ø
x Î[Ãlv, Ãsl ].
Hs = cpsTs
ps
=
uc2ρ0s
èççççæèçççæ
ρ ρ0s
-1ø÷÷÷ö÷
+
βs
(T
-T0
)÷ø÷÷÷ö
x Î[Ãsl, xs ].
(4)
(5) (6)
В рассматриваемой системе уравнений
использованы следующие обозначения и единицы: ρ [г/см3] – плотность, u [см/с],
р [бар] – гидродинамическая скорость и давление, uс[см/с] – скорость звука, H [Дж/см3] – энтальпия, W [Вт/см2] – тепловой поток,
α [1/см] – коэффициент поглощения лазерного излучения, G [Вт/см2] – плотность мощно-
сти лазерного излучения, cp [Дж/г K] – удельная теплоемкость при постоянном давлении,
β [1/K] – температурный коэффициент линей-
ного расширения, T0 [K] – комнатная темпе ратура, Tm [K] – равновесная температура плавления, ρ0, s, ρ0, l – начальные плотности твердой и жидкой фазы, задаваемые соответ-
ственно при температурах T0 и Tm. Индексы sur, s, l, υ обозначают принадлежность ве-
личин соответственно к поверхности, твердой,
жидкой и парообразной средам, а индексы sl
и lυ – принадлежность к межфазным грани-
цам раздела твердое тело–жидкость x = Γsl и границе испарения Гlυ. Система уравнений
80 “Оптический журнал”, 78, 8, 2011
(1)–(6) дополнялась граничными и начальными
условиями.
Левая граница х = х0, полагалась закрепленной и теплоизолированной (рис. 1).
x = x0:
λ
¶Ts ¶x
= 0,
us
= 0.
(7)
На межфазной границе x = Γsl формулируется модель поверхностного плавления, состоя-
щая из трех законов сохранения, дополненных
условием непрерывности температуры фазово-
го перехода.
x = Ãsl: ρsusl = ρl (us -ul + usl ), ρsu2sl + ps = ρl (us -ul + usl )2 + pl,
(8) (9)
λl
¶Tl ¶x
- λs
¶Ts ¶x
=
Lnme ρs usl ,
(10)
Ts = Tl = Tsl = Tm,
(11)
где
( )( )Lnme = Lm + cpl - cps
Tsl -Tm, 0
+
ρs ρs
+ -
ρl ρl
(us
- ul )2
2
– неравновесная теплота плавления, Lm [Дж/г] – равновесная теплота плавления, λ [Вт/см K] –
коэффициент теплопроводности. На правой облучаемой границе x = Γlυ(t) фор-
мулируется модель поверхностного испарения,
записанная в приближении кнудсеновского
слоя и состоящая из трех законов сохранения
и двух дополнительных условий, характери
зующих степень неравновесности фазового пе-
рехода:
x = Ãlu : ρlulu = ρu (ulu + ul -uu ),
(12)
pl + ρlu2lu = pu + ρu (ulu + ul -uu )2 ,
(13)
WlT = λl
¶Tl ¶x
= ρluluLnue
+ σTs4ur ,
(14)
Tu = TlαΤ (M), ρu = ρsatαρ (M),
(15)
Glu
=
A(Tsur
)G0
expèççççæç-4çèçæç
t tL
ø÷÷÷÷ö2
÷ø÷÷÷÷÷ö,
(16)
где
Lnue
=
Lu
(Tl
)+
cpu
(Tl
-Tu)+
ρl ρl
+ -
ρu ρu
(ul
- uu)2
2
– неравновесная теплота испарения, Glυ, G0 – поглощенная поверхностью энергия лазер-
ного излучения и ее максимальное значение,
Asur [%] – поглощательная способность поверхности, τL [нс] – полуширина лазерного импульса, Lυ [Дж/г] – теплота кипения, M = = ulυ/uc – число Маха, uc = (γRTυ)1/2 – скорость звука на внешней стороне кнудсеновского слоя,
γ = cp/cυ, R [Дж/г K] – универсальная газовая постоянная, αT(M), αρ(M) – коэффициенты Крута [7, 8], при М = 1 αT(M) = 0,633, αρ(M) = 0,326.
Давление насыщенного пара вычисляли, ис-
пользуя уравнения Клапейрона–Клаузиуса:
psat
=
pb
expæççççè
Luµ RTb
æçççè1-
Tb Tl
ö÷÷÷ø÷÷ö÷÷÷ø,
ρsat = psat /RTl,
(17)
где pb, Tb – соответственно давление и температура кипения при нормальных условиях. В на-
чальный момент времени t0 температуру и плотность полагали постоянными, а скорость,
равной нулю:
Ts = T0, ρs = ρ0,s, us = 0.
(18)
Оптические и теплофизические характеристики алюминия и кремния
На форму оптоакустических сигналов большое влияние оказывают теплофизические и оптические характеристики материала мишени [9, 10]. Используемые в расчетах характеристики материалов, представленные на рис. 2a–б и 3, взяты из справочных данных [11, 12]. Вертикальными линиями на рисунках отмечены равновесные температуры плавления и испарения каждого из материалов. Все теплофизические и оптические характеристики обоих материалов претерпевают разрыв при переходе через значение равновесной температуры плавления (на рис. 2, 3 она отмечена вертикальными линиями). Учитывая, что один из рассматриваемых материалов является металлом, а другой полупроводником, их теплофизические и оптические характеристики сильно различаются. При переходе через фронт плавления теплофизические характеристики алюминия скачком уменьшаются, а кремния – возрастают. В таблице приведены значения теплофи зических параметров, характеризующие свойства рассматриваемых материалов.
Алгоритм решения. Основная сложность решения задачи Стефана заключается в наличии двух подвижных границ, положение которых неизвестно и должно определяться в ходе расчетов, при этом размеры твердой и жидкой подобластей в процессе расчетов могут меняться на несколько порядков. Анализ вклада
“Оптический журнал”, 78, 8, 2011
81
, Вт/см3 , г/см K
3
2
1
Tm
1000
, Вт/см3 , г/см K
(а)
1
Tu
3000
5000
(б)
2 3
7000
Ср, Дж/гK
3 2 1
Ср, Дж/гK
31
3
22
3 2
11
1000
Tm
2000
Tu 3000 4000 Т, K
Рис. 2. Теплофизические характеристики алюминия (а) и кремния (б). 1 – плотность; 2 – те-
плоемкость; 3 – теплопроводность; Тm – температура плавления; Тν – температура испа рения.
фазового перехода в оптоакустический сигнал требует точного определения положения межфазной границы и ее скорости. Наиболее эффективным методом для решения поставленной задачи является метод динамической адаптации [13–16], позволяющий производить расчет с явным выделением межфазных границ.
Метод динамической адаптации. Численное решение задачи (1)–(18) осуществляется методом динамической адаптации, в основу которого положена идея перехода к произвольной нестационарной системе координат, осуществляемого с помощью искомого решения. В про-
100
A, %
80
60
(а)
40
20
0
Tm
1000
Tu
3000
5000
7000
100
A, %
80
60
(б)
40
20
Tm
0
1000
2000
3000 Т, K
Рис. 3. Поглощательная способность алюминия (а) и кремния (б). Тm – температура плавления. Тv – температура испарения.
извольной нестационарной системе координат проблема описывается расширенной дифференциальной системой уравнений, часть из которых описывает физическое явление, а другая часть – динамику узлов расчетной сетки. Проблемы, связанные с подвижными границами, снимаются путем перехода к произвольной нестационарной системе координат, в которой узлы сетки и границы оказываются неподвижными. Преобразование координат осуществляется автоматически с помощью искомого решения, что позволяет производить размещение узлов сетки в зависимости от особенностей решения.
Переход из физического пространства с переменными (x, t) в расчетное с произвольной нестационарной системой координат (q, τ) осу-
Теплофизические параметры алюминия и кремния
Минерал
Атомный вес, г/моль
Al 26,98 Si 28
Температура
Tm K
933 1683
Tb K
2793 3514
Коэффициент теплового расширения, 1/K
Твердая фаза βs
2,33×10–5 4,65×10–3
Жидкая фаза βl
5,00×10–5 1,40×10–4
Теплота перехода
Плавление Lm (Дж/г)
400,30 1797,0
Кипение Lb (Дж/г)
10860 13720
82 “Оптический журнал”, 78, 8, 2011
ществляется с помощью замены переменных общего вида
x = f(q, t), t = t,
(19)
имеющей обратное невырожденное преобразо-
вание q = j(x, t), t = t. Частные производные в
системе координат (q, τ) имеют следующий вид:
¶ ¶x
=
1 ψ
¶ ¶q
,
¶ ¶t
=
¶ ¶t
+
Q ψ
¶ ¶q
,
(20)
где ∂x/∂τ = –Q – скорость движения новой си-
стемы координат относительно исходной, под-
лежащая в дальнейшем определению.
В нестационарной системе координат две
ф азовые подобласти [Гlν, Гsl] ∪ [Гsl, xs] с по движными границами Гlν, Гsl отображаются на две расчетные подобласти [0, qsl] и [qsl, qs] с неподвижными границами. До появления
жидкой фазы физическая область [Гlν, xs] отображается на расчетную подобласть [0, qs]. Физическая координата х в расчетном простран-
стве становится новой неизвестной функцией.
В новой системе координат динамика расчет-
ной сетки описывается дополнительным диф-
ференциальным уравнением:
êëéê
¶ψ ¶t
=
-
¶Q ¶q
úúùû k
k = s, l,
(21)
где ψ = ∂x/∂q – метрический коэффициент, Q – функция, определяющая конкретный вид
преобразования в соответствии с особенностями
задачи.
Используя преобразование переменных (19),
запишем дифференциальную задачу (1)–(17)
в произвольной нестационарной системе координат (q, τ):
êëêé
¶(ψρ)
¶t
+
¶ ¶q
(ρ(u
+
Q))=
0ùûúú
,
k
êëêé
¶ ¶t
(ψρu)+
¶ ¶q
(ρu(u
+
Q))=
-
¶p ¶q
ùúúû
,
k
éêêêë
¶ ¶t
æçççèçψρæççççèH
+
u2 2
ö÷ø÷÷÷ö÷ø÷÷÷÷
+
¶ ¶q
æçççèçρ(u
+
Q)æççççè
H
+
u2 2
÷ö÷÷÷ø÷÷÷÷÷öø
=
=
-
¶ ¶q
(
pu)+
¶ ¶q
ççèçæç
λ(T)
ψ
¶T ¶q
÷ø÷ö÷÷ûúúúù
,
k
(22)
êéëê
¶G ¶q
+
αψG
=
0úùúû
,
k
q Î[0, qsl ][qsl, qs ], t Î[t0, tend ], k = s, l.
Граничные условия. q = q0: u = W = 0, Q = 0;
(23)
q = Γsl: ρs (Qsl + us )= ρl (Qsl + ul ), ps + ρs (Qsl + us )2 = pl + ρl (Qsl + ul )2,
(24) (25)
Wl - Ws = ρsuslLnme, usl = Qsl + us,
q = Γlu: ρl (Qlu + ul )= ρu (Qlu + uu ), pl + ρl (Qlu + ul )2 = pu + ρu (Qlu + uu )2,
(26) (27) (28)
WlT
=
WlT
-
λl ψl
¶Tl ¶q
= ρluluLnue + σTs4ur ,
ulu = -(Qlu + ul ),
Glu
=
A(Tsur
)G0expçèçççæç-4èçççæ
t tL
÷÷ø÷÷ö2
÷÷ø÷÷ö÷÷.
(29)
t = t0: Ts = T0, us = 0, ρs = ρ0, s, ψ = 1. (30)
Таким образом, при переходе к произвольной нестационарной системе координат исходная дифференциальная модель трансформируется в расширенную дифференциальную систему (21)–(30), в которой появляется дополнительное уравнение типа (21), являющееся уравнением обратного преобразования. Его тип, свойства и вид краевых условий зависят от конкретного вида функции Q. Для построения равномерных (квазиравномерных) на каждый момент времени сеток в областях с подвижными границами функция Q задается в виде [13, 14]:
êéêëQ
=
-D
¶ψ ¶q
ûùúú
,
k
k = s, l,
где коэффициент диффузии D выражается че-
рез геометрические и скоростные параметры
задачи: D = L2(t)(|υsl| + |υlυ|)/Δx. Разностная аппроксимация. Численное ре-
шение дифференциальной модели (21)–(30)
осуществлялось при помощи конечно-разност-
ного метода, согласно которому уравнения ги-
дродинамики аппроксимировали семейством
консервативных разностных схем, получен-
ных интегро-интерполяционным методом [17].
Для построения семейства разностных схем в расчетном пространстве вводится сетка wqt с неравномерными шагами hk, i, τj по пространственной q и временной τ переменным:
ω = {ωl ∪ ωs}×{ωτ},
“Оптический журнал”, 78, 8, 2011
83
где
wl
=
îïïìïíïiq=l,
i, ql, i+1/2; 0,.., Nl -1,
ql, i+1 = ql, ql, 0 = 0,
i+ qNl
hl, =
i+1, qsl ,
ql, i+1/2 hl, 0 = 0,
= ql, i hl, Nl
+
+1
0,5hl, =0
i+1
ïüïþïýï,
ws
=
ìïïïîïíiq=s,
i, 0,
qs, ..,
i+1/2; Ns -1,
qs, i+1 = qs, i + hs, i+1, qs, 0 = 0, qNs = qsl,
qs, i+1/2 hs, 0 = 0,
= qs, i hs, Ns
+
+1
0,5hs, =0
i+1
ýþüïïïï,
{ }wt = tj, tj+1 = tj + Δtj+1, j = 0,..., J, t0 = t0 .
При разностной аппроксимации уравнений
к узлам с целочисленными индексами отно
сятся сеточные функции p, Q, x, к точкам с полуцелыми индексами (qm+1/2, τj) – сеточные функции ψ, ρ, u, T. Аппроксимация системы
уравнений (22) производилась с помощью полностью неявной разностной схемы, имеющей порядок аппроксимации Q(τ + h).
Для дифференциальной системы (21), (22) семейство разностных схем имеет вид:
ïíïîïïïïì(ψρ)ij++11/Δ2t-j+(1ψρ)ij+1/2
+
ëéρ(u
+
Q)ûùij++11 - ëéρ(u
hi+1
+
Q)ûù
j+1 i
=
0ýïïïïïþüïk,
(31)
ïíîïïïïïì(ψρu)jj++11/Δ2t-j+(1ψρu)jj+1/2
+
ëéρu(u
+
Q)ûù
j+1 i+1
-
ëéρu(u
+
Q)ûù
j+1 i
hi+1
=-
pij++11 + pij+1 hi+1
ïþïïïýïïük,
(32)
( ) ( ) ( ) ( )ïíïïïïïïïïîïïïïïïïì=éêëψ-ρ
ε + u2
2
ùúû
j+1 i+1/2
-
éëêψρ
ε
Δtj+1
(pu)ij++11 +(pu)ij+1 + 1
hi+1
hi
+ u2
2
ùûú
j i+1/2
+ éêêëρ(u + Q)
ε + u2
2
ij++11ùúúû - éêëêρ(u
hi+1
( )ççççèçæ
λij++11
ψ
j+1 i+1
Tij++31/2 -Tij++11/2 0,5 hi+1 + hi
-
λij+1 ψij+1
T0i,j+5+1(1/h2i-+Thi-ji+-111/2)ø÷÷÷÷÷ö
+
Q)
ε + u2
2
ij+1ùúúû
=ïïýïïïïïïïïïïüïþïïïk,
(33)
éêëêê
ψij++11/2 - ψij+1/2 Δtj+1
=
-
Qij++11 - Qij+1 hi+1
ùúûúú
,
k
(34)
ëêêêé
Gij++11 - Gij+1 hi+1
+
α(ψG
)j+1
i+1/2
=
0úúûúù
,
k
(35)
éêêêëQij+1
=
-Dij
ψ
j+1 i+1/2
-
ψij-+11/2
(hi + hi+1 ) 2
úùúûú,
ëêêêé
xij++11 - xij+1 hi+1
=
ψij++11/2
ûúùúú
,
k
где i = 0, …, Nl – 1 при k = l, i = 0,…,Ns – 1 при k = s.
Аппроксимация граничных условий (23)–
(30) выполнялась следующим образом:
84
q = q0:
ψsj,+01 λsj,+01
Tsj,+01 -Tsj,+1/12 hs, 0 /2
= 0,
us, 0 = 0, Qs0 = 0,
(36)
q = Γsl : (ρsVsl )j+1 = ρlj+1 (us - ul + vsl )j+1, (37)
( )ρsVs2l j+1+ psj+1=
( )= ρlj+1
(us -ul + vsl )j+1
2
+
plj+1.
(38)
λs, N-1 ψs, N-1
Tsj,+N1-1/2 -Tsj,+N1-1 hs, N-1/2
-
λl, 0 ψl, 0
Tl,j+01 -Tl,j1+/12 hl,1/2
=
= Lm (ρsvsl )j+1,
(39)
Ts = Tl = Tm.
“Оптический журнал”, 78, 8, 2011
q = Γlu:
(ρvvlv )j+1= (ρv (ulv -uv + vlv ))j+1,
( )Qljv+1 = - vljv+1 + uljv+1 ,
(40)
( )ρlvvl2v j+1+ pljv+1=
( )= ρv (ulv -uv + vlv )2 j+1+ pvj+1,
(41)
λsj+ur1
( )ψsj+ur1
Tkj,+11/2 -Tsju+r1 hk,1/2
= (ρsur LvVkv )j+1+
σT4
j+1,
(42)
( )(Glu )j+1= A Tsju+r1
G0
expçèççæçç-4çèæçç
t tL
÷ø÷÷÷ö2
÷÷ø÷÷÷÷öj+1.
(43)
Величины в полуцелых узлах для сеточных функций, отнесенных к целым точкам, вычисляются как полусумма значений в ближайших узлах. При определении величин в узлах с целочисленными индексами для функций, о тнесенных к полуцелым точкам, используется линейная интерполяция:
ym
=
ym-1/2ψm+1/2 ψm+1/2
+ ym+1/2ψm-1/2 + ψm-1/2
.
(44)
Результаты моделирования
Рассмотрим и проанализируем изменения амплитудных характеристик оптоакустических сигналов в алюминии и кремнии, вызванных короткоимпульсным лазерным воздействием.
Моделирование формы акустического си гнала проводилось для лазерного импульса гауссовой формы с длительностью на полувысоте τL = 3 нс, в диапазоне пиковых интенсивностей G0 = 107 – 1010 Вт/см2. Кривые давления получены для глубины хℓ =10 мкм. Точка расположена вне зоны теплового влияния импульса и на рассматриваемом интервале времени основной сигнал в ней не искажается отраженным. Выбранный диапазон интенсивности излучения для данной длительности τL соответствует области доминирования процесса плавления обоих материалов, а испарение еще не вносит существенного вклада. Сравнительный анализ оптоакустических сигналов должен дать представление:
–– о пороге плавления обоих материалов; –– о характерном профиле оптоакустических сигналов для металлов и полупроводников; –– о влиянии теплофизических и оптических параметров облучаемых материалов на динами-
ку фазовых переходов и амплитуды оптоакустических сигналов.
Для теплового механизма генерации звука существуют оценки [1], полученные из совместного решения линейных уравнений теплопроводности и гидродинамики, согласно которым:
P(t)=
χβρ0
¶Tsur ¶t
+
1 cpα
¶G¶t(t),
χ
=
λ cpρ0
.
(45)
На рис. 4а–в представлены временные за-
висимости температуры поверхности, скоро-
сти фронта плавления и давления для крем-
ния. Моделирование показало, что, несмотря
на более высокую, чем у алюминия, темпе
ратуру плавления TmSi (1683 K)> TmAl(933,6 K),
порог плавления кремния оказался ниже
(рис. 4a). Пороговая интенсивность лазерного воздействия составила G0 = 2,4×107 Вт/см2. Температура плавления на облучаемой поверх-
ности кремния достигается на заднем фронте
лазерного импульса (рис. 4a). Быстрому нагре-
ву приповерхностных слоев кремния способству-
ет значительное уменьшение теплопроводности
в твердой фазе с ростом температуры. Появле-
ние жидкой фазы характеризуется скоростью Vsl (рис. 4б), отрицательная ветвь которой характеризует плавление, а положительная – за-
твердевание. Оптоакустический сигнал явля-
ется динамической характеристикой, которая
весьма чувствительна к темпу ввода энергии в
систему. На начальной стадии нагрева импульс
давления (рис. 4в), положителен. Величина его
растет с ростом температуры до значения рa ≈ 3,5 бар, а затем замедляется из-за уменьшения производных ∂Tsur/∂t, ∂G(t)/∂t и плавно спадает после смены их знака. Момент обра-
зования жидкой фазы характеризуется отрицательным скачком плотности Δρ = ρs – ρl 0, отрицательная полуволна – сжатию ∂Ts/∂t 0. Это проявляется на развертке давления в виде резкого коротко-
го скачка вверх, вплоть до нулевого значения давления рa. Последующее затвердевание происходит при быстро уменьшающемся энергов-
кладе, в результате чего максимальная скорость
затвердевания (≈ 60 м/с) в 1,5 раза превышает
максимальную скорость плавления (рис. 5б), и в
7,5 раз превышает максимальную скорость за-
твердевания кремния. Столь высокая скорость
затвердевания вызывает резкое увеличение отрицательного давления, рa ≈ 15 бар (рис. 5в).
При повышении интенсивности излучения до G0 = 2×108 Вт/см2 оптоакустические сигна-
86 “Оптический журнал”, 78, 8, 2011
Тsur, K
800
400
0 – 4 – 2
0
(а) usl, м/с
(б) Р, бар
40 10
20 0
0
– 20 – 10
– 40 – 20
2 4 6 t, нс
– 4 – 3 – 2 –1 0 1 2 3 t, нс
– 4 – 2 0
2
(в) 4 6 t, нс
Рис. 5. Зависимость от времени температуры поверхности Тsur, скорости плавления usl и давления р для алюминия. Пунктиром обозначена зависимость поглощенной энергии лазерного излучения от времени. Полуширина лазерного импульса 3 нс. G0 = 3,5×107 Вт/см2.
Тsur, K
2500 2000
1500
(а)
usl, м/с
300
200
100
0
1000
– 100
500 – 200
0 – 2
0
2
4
– 300
6 t, нс
– 4
04
(б) ulv, м/с
8 10–3 4 10–3 0
– 4 10–3 – 8 10–3 – 1,2 10–2
8 t, нс
Р, бар
200 100
0 – 100
– 2
0
(в) 2 4 t, нс
Рис. 6. Зависимость от времени температуры поверхности Тsur, скорости плавления usl, испарения ulν и давления р для алюминия. Пунктиром обозначена зависимость поглощенной энергии от времени. Полуширина лазерного импульса 3 нс. G0 = 2×108 Вт/см2.
Тsur, K
2400
1600
800
0 – 4 – 2
0
(а)
24
usl, м/с
2 10–3
(б)
0
– 2 10–3
– 4 10–3
– 6 10–3
– 8 10–3
– 1 10–2
6 t, нс
– 4 – 2 0 2 4
ulv, м/с
4 10–5
Р,
бар
10
0 – 4 10–5
0
– 8 10–5
– 10
– 1,2 10–4
– 20
– 1,6 10–4 – 2 10–4
6 t, нс
– 30 – 40
– 4
– 2
0
2
(в) 4 t, нс
Рис. 7. Зависимость от времени температуры поверхности Тsur, скоростей плавления usl, испарения ulν и давления р для кремния. Пунктиром обозначена зависимость поглощенной энергии от времени. Полуширина лазерного импульса 3 нс. G0 = 2×108 Вт/см2.
лы алюминия и кремния качественно меняются (рис. 6a–в, 7а–в). Момент плавления перемещается на передний фронт лазерного импульса, где возрастает роль поглощательной способности материала. В алюминии поглощательная способность при плавлении скачком увеличивается, усиливая скорость плавления до 300 м/с (рис. 6б) и, как следствие, происходит скачок давления до 270 бар (рис. 6в). Поверхность при этом разогревается до температуры TsAulr ³ TbAl и становится заметным вклад испарения, скорость которого достигает vlν ≈ 8×10–3 м/с
(рис. 6б). Амплитуда оптоакустического импульса после всплеска давления, вызванного скачком поглощательной способности и положительным скачком плотности, под совместным воздействием плавления, испарения и уменьшения интенсивности излучения быстро спадает, переходя в отрицательную полуволну.
Оптоакустический сигнал в кремнии имеет еще более сложную временную структуру (рис. 7а–в). Для кремния поглощательная способность поверхности скачком уменьшается
“Оптический журнал”, 78, 8, 2011
87
(рис. 7а), что в совокупности с отрицательным скачком плотности приводит к увеличению отрицательного скачка давления до 35 бар (рис. 7в). Продолжающийся нагрев позволяет уменьшить отрицательное давление до нуля, которое затем снова увеличивается до 10 бар, из-за совместного воздействия испарения, скорость которого (vlv ≈ 4,5×10–3 м/с, рис. 7б) становится заметной, и уменьшения интенсивности лазерного импульса. Процесс затвердевания приводит к формированию положительной полуволны с максимальным значением pa ≈ 10 бар. Из-за более низкой п оглощательной способности темп нагрева кремния ниже, чем у алюминия, чем и объясняются более низкие максимальные значения температуры поверхности, скорости плавления и скачков давления.
Заключение
Выполненный с помощью математического моделирования сравнительный анализ опто акустических сигналов, вызванных лазерным воздействием на алюминий и кремний, показал, что:
1. Плавление существенно меняет вид оптоакустического сигнала в кремнии и алюминии.
2. Скачок давления, связанный с плавлением, положителен у алюминия (Δp > 0) и отрицателен у кремния (Δp