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

УДК517.11+519.62+536.21(45)

ТИТОВ Александр Константинович, кандидат физико-математических наук, доцент кафедры информатики, вычислительной техники и методики преподавания информатики Поморского государственного университета имени М.В. Ломоносова. Автор более 20 научных публикаций

ЕШЕВСКИЙ Олег Юрьевич, кандидат физико-математических наук, доцент кафедры физики Поморского государственного университета имени М.В. Ломоносова. Автор 40 научных публикаций

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

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

Математическое моделирование, численное моделирование, компьютерное моделирование, лазерное испарение

Моделирование действия лазерного излучения (ЛИ) на дисперсные частицы является актуальной задачей как для многочисленных технических приложений, так и в области фундаментальных исследований. Воздействие лазерного излучения на среды, в частности на металлические частицы, сопровождается различными процессами, включая поглощение, нагрев, плавление, испарение, разлет разогреваемого излучением пара, образование плазмы и т.д. Ранее одним из авторов [1,2] рассматривалось действие ЛИ на паровой ореол мелких металлических частиц в стационарном режиме. Процессы распространения тепла и плавление частицы в указанных работах не рассматривались. В работах [3, 4] предложены модели плавления металлической плоской преграды на основе задачи Стефана. В отличие от [4], в [3] предлагается энтальпийная формулировка, позволяющая не рассматривать непосредственно

движение границы раздела фаз. В предлагаемой работе будем рассматривать распространение фронта плавления по сферической частице в соответствии с [4].

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

Поглощение лазерного излучения. Как показано в [5], значительная часть излучения в ИК-диапазоне, падающего на металлическую частицу, будет отражаться и только К - 0,01 этого излучения поглощается частицей. При этом поглощение будет происходить в тонком приповерхностном слое (скин-слое). Толщина скин-слоя для излучения С02 - лазера (длина волны X = 10,6 мкм) в алюминии д - 1 мкм. Будем моделировать действие лазерного излучения интенсивностью /на частицу радиальным распределением источников тепла в шаровом слое частицы радиусом го толщиной д с плотностью

/{г) - • х(г - г0 + 8), (1)

о

где х(х) ~ единичная функция Хевисайда.

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

/(г)

2 1 д 2 ди ~ ч ди 2 К ак ——г — +/к (г) =—, ак = г дг дг д1 ск ■рк

к , Л (Г) =-

ск 'Рк

-к,

ди

' дг

г=т ,+0 рп

ди +щ — дг

=м-р-К

рк

г^г ,-0 рп

к=1,2.

(2)

Здесь кк,ск, рк и ^ - коэффициент теплопроводности, удельная теплоемкость, массовая плотность и удельная теплота плавления, соответственно; индекс к= 1 относится к твердой фазе, а к=2 - к жидкой, индекс рк - к границе раздела фаз; ¥рк - скорость движения поверхности раздела фаз.

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

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

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

мент времени она одинакова и равна нулю всюду в частице, т. е. и(г, 0)=0.

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

д и д г

= 0.

(3)

Разумно также предположить, что в центре частицы тепловой поток также отсутствует:

д и д г

= 0.

(4)

г-0

Однако для последнего условия, как будет показано ниже, возникают проблемы аппроксимации разностным уравнением. Рассмотрим эту проблему позднее.

На границе раздела фаз температура равна температуре плавления Тръ.

«(гРи,{) = иРи = Тръ - Т^= сош1- (5)

Математическая модель. Соберем все уравнения, рассмотренные выше (1)-(5), и окончательно сформулируем математическую модель.

Первый этап - до момента Ь начала плавления

3 и

Ъ ~Т^Гг2^Г + Мг) = >0 < г < г0’° < t < ?1>

о t

2 1 3 2 & и

~2----Г '

г д г д г

д и 8 г

= 0, и{г,0) = 0.

Второй этап - плавление.

а1

1 д 7 ди

ди

г +Кг)=-г. 0 <г < гРк;

г дг дг

і4гф ,і) = иф = сопії, -к.

+К

дґ ди 1 дг

рк

рк

2 1 5 , Зм ~ 5м «г —■— г — +/2(г) =—, <г<г0, ^ <? <?а;

г ог дг от

і)=9<г),0 < г < г0,

=О,

ди

дг

=0.

(7)

Здесь ф(г) - результаты счета первого этапа для момента времени

Следует добавить к системе (7) уравнение движения границы раздела фаз:

рк

(і і

V

рк *

(8)

Математическая модель (6)-(8) позволит:

1) проследить динамику температуры частицы и плавления частицы;

2) определить времена нагрева до температуры плавления и полного плавления частицы;

3) определить скорость движения фронта плавления в различные моменты времени.

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

^ и г ^

и ----, х - —, т - —.

Г, Г0 іп

' рк

(9)

Масштаб времени іо выберем, положив в уравнении теплопроводности аі2іо/го2=1:

го • С1 • Л К,

Тогда задача распространения тепла на первом этапе примет вид:

1 д 2 ^ и х2 д х д х д и д х

= о,

д х

дт

= 0, и (х,0) = 0.

Второй этап:

1 д і д и

(10)

д и

——х — + С ■ Дх) = —, 0 < х < хрк; X О X С X

и (хрк ,т) = 1, -К.

ди д х

+ к,

дт д и 1 д х

'рк

дх . рк &

я, 1 д і д и

д и

Х Т- + СХ • ДХ) = Т“ > ХрИ < Х < !> ;

а1 х д х д х

и (х,Т[) =ф(х),0 < х < 1,

д и д х

дт'

= о,

д и д г

= 0

В уравнениях (10)-(11):

(11)

С = -

*х -8- {ТрН - )

к\ ‘ (0 ‘ (Трк — )

/(^) = х(х -1 + —), С1 = С

Рх ■

р2 • С2

Рх • И- гс

=■

УрН ‘ ¿0

(12)

Для численного моделирования задачи (10)-(12) использовалась явные конечно-разностные схемы.

Конечно-разностная схема. Уравнения теплопроводности в (10)-(11) для однофазного режима или твердой фазы двухфазного режима и внутренних узлов сетки аппроксимировались консервативной [6] явной конечно-разностной схемой:

{х, + к/2)2 {х, - к/2)2

+ С • /, . (13)

Для жидкой фазы:

+1, і - , І

1 2 2 п • а,

(X,. + к /2)2 (X,. - к ну

+Су/,.

(14)

Разностная схема для краевого условия второго рода на поверхности частицы может быть получена путем добавления дополнительного N+1 узла по отношению к N узлу поверхности [7]. В итоге краевое условие для х^=1 в (10) аппроксимируется разностной схемой:

У ]+1, N V], N

(хк + к И)2 (хК - к/2)

V V-'- N

+ С ■ /• •

(15)

Аналогично аппроксимируется и условие для (11).

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

л 2г2 ^ ЗгЗ ^ ^

] ■ 4-п-г0 • к = с, ■ рх ---п-г0 • к

З о Ї

Или

к,

V! ,1 - V! ,0

= С1 • Р1 • — го ■ к ■

Окончательно имеем:

^,1 - V,,0 1 V,-+1>0 - V,->0

(16)

На шаг временной сетки т накладывается ограничение в виде условия Куранта [6, 7]:

т < ■

к:

которое является достаточно жестким и существенно увеличивает объем и время вычислений. В таких случаях [6, 7] принято использовать неявные или явно-неявные схемы с весами, например схему Крэнка-Николсона. Тем не менее в среде МаШСаё была реализована компьютерная модель по алгоритмам (13)—(16), которая позволила выполнить вычисления и получить результаты.

Результаты вычислений и анализ. Ниже на рисунке приведены результаты вычислений для частицы алюминия радиусом 10 мкм в поле излучения С02-лазера интенсивностью 108 Вт/см2. На рисунке горизонтальной пунктирной линией отмечена температура плавления (7^=933 К). Температурные профили {То - для момента 2,0-10-8 с, Тг - 6,0-10-8 с, Т2 - 8,6-10-8 с, Т3 - 2,2-10"7 с, Т4 - 3,4-10"7 с, Т5 - 6,1-10'7 с ) показывают динамику нагрева частицы. По достижению температуры плавления в приповерхностном слое появляется поверхность раздела фаз, которая движется в глубь частицы. Как видно из графиков, перед фронтом волны плавления температура выравнивается быстрее, чем движется сам фронт. Практически она не меняется от поверхности раздела фаз и до центра частицы (равна температуре плавления). Частица, пока не расплавилась полностью, представляет сбой твердое ядро, окруженное расплавленным металлом. Средняя скорость движения фронта плавления, как мы видим, порядка 15 м/с. К концу движения она достигает -70 м/с. Можно выполнить оценку полного плавления частицы, сложив энергозатраты на нагрев до температуры плавления всей частицы, на плавление, нагрев жидкой фазы, и затем поделить на поглощаемую мощность. Такая оценка дает почти совпадающий с полученным значением (6,1-Ю7 с) результат, что свидетельствует об адекватности модели. Полученные результаты показывают, что при полном плавлении данной частицы ее жидкая фаза прогре-

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

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

Н

То

Ti

Т2

Т3

Т4

т,

Радиальная координата, м

Динамика нагрева и плавления частицы

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

1. Интенсивное испарение коллектива частиц с учетом тепловой и ионизационной неравновесностей под действием излучения / А.Г. Лесскис, А.К. Титов, A.A. Юшканов, Ю.И. Япамов II Теплофизика высоких температур. 1991.Т. 29,№5. С. 864-871.

2. Лесскис А.Г., Титов А.К., Юшканов A.A. Установившееся радиальное движение разогреваемого излучением пара в вакуум от интенсивно испаряющейся частицы металла II Теплофизика высоких температур. 1995. Т. 33, № 4. С. 578-582.

3. Моисеев К.В., Кузенов В.В. Численный анализ одномерной задачи Стефана при лазерном воздействии на металлические преграды II Физико-химическая кинетика в газовой динамике. URL: www.chemphys.edu.ru/pdf/2008-09-01-042.pdf.

4. Суровцев А.Н., Суровцева H.A., Титов А.К. Нагревание и плавление металлической поверхности лазерным излучением// Физ. вести. Помор, ун-та: сб. науч. тр. Вып. 7. Архангельск, 2008. С. 42-46.

5. Хюлст Г. ван де. Рассеяние света малыми частицами. М., 1961.

6. Самарский А.А. Теория разностных схем. М., 1989.

7. МарчукГ.И. Методы вычислительной математики. М., 1989.

Titov Aleksandr, Eshevsky Oleg

NUMERICAL MODELING OF HEATING AND MELTING OF SPHERICAL METAL PARTICLE BY LASER EMISSION

The mathematical model of metal particle’s heating and melting process by laser emission based on Stefan’s challenge is given. The numerical scheme of equating is developed. Evaluation and analysis of results are carried out.

Контактная информация: Титов Александр Константинович e-mail\ titov.aleksandr@pomorsu.ru Ешевский Олег Юрьевич e-maiI: eshevsky.oleg@pomorsu.ru

Рецензент- Попов В.Н., доктор физико-математических наук, доцент кафедры общей физики Поморского государственного университета имени М.В. Ломоносова