УДК 621.375
В. И. Букатый, А. А. Попов, А.М. Шайдук
Компьютерное моделирование процесса горения и испарения углеродной частицы в поле лазерного излучения
В связи с развитием теоретических и практических исследований в области взаимодействия мощного лазерного излучения (МЛИ) с веществом, а также компьютерной техники возникла необходимость численного моделирования процессов, проходящих в веществе, в частности, в дисперсных системах при воздействии лазерного излучения. В настоящее время существует ряд работ, например [ 1, 2 ], посвященных воздействию МЛИ на углеродные частицы, однако отсутствие в достаточной мере экспериментальных данных и соответствующих теоретических расчетов не позволяет адекватно построить модель взаимодействия.
Целью данной работы является компьютерное моделирование процессов, проходящих в углеродных частицах микронных размеров под действием МЛИ в широком диапазоне интенсивностей излучения.
Задача заключается в следующем: есть система дифференциальных уравнений, описывающая поведение частицы в поле МЛИ, т. е. изменение ее радиуса и температуры с течением времени под действием поля МЛИ, необходимо написать компьютерную программу, которая бы позволяла сравнительно просто получать характеристики процесса (параметры частицы и окружающей среды!) с заданными начальными условиями. Уточним задачу. Пусть имеется сферическая твердая частица и на нее постоянно падает заданной мощности излучение. Частица задается начальным радиусом, температурой и константами, характеризующими ее как вещество, такими как теплота испарения, плотность, масса и др. Необходимо определить, что будет с этой частицей, т.е. каков должен быть радиус частицы, если на нее воздействовать излучением заданной интенсивности, такой, чтобы температура ее поверхности была определенной, как долго будет поддерживаться данная температура и в течение какого времени? Система дифференциал ьных уравнений состоит из двух уравнений, которые и описывают поведение частицы в поле МЛИ.
da _ _ Ks (a,Ts) dt p
m
4pa 2p J3'
dT 3 )
— _----------3-------F(a, T),
dt 4pa pCp
(1)
Начальные условия: a ( 0 ) _ a 0 ; Ts ( 0 ) _ T0 . В систему уравнений (1), (2) входят следующие выражения:
F (a, Ts ) _ Я„огл + Чгор ~ Ч тепл Уизл Чисп ’
m
Ks (a, Ts) _ 4p 2 p
4pa p
Чпогл _Pa2Kn (a> 1) 1;
Vi,
Чгор = 4ра 2 <2эфф (а, Т )Ks (а, Ts );
Тз
Я тепл = 4жи\Х(Т) йТ;
То
Яшш = Ала2а(Т? - Т04);
Чисп = ^03 ,
где ^, ^ - поток кислорода к поверхности частицы и поток вещества от частицы соответственно; у -эффективный коэффициент стехиометрии комплекса реакций р ; Ср - плотность и теплоемкость вещества; т - масса вещества; Кп - фактор эффективности поглощения излучения частицей; Сэфф
- эффективный тепловой эффект комплекса гетерогенных реакций горения; 1 - теплопроводность воздуха; и - постоянная Стефана Больцмана; Т0
- температура среды на бесконечности; Тз - температура поверхности частицы; Ь - теплота испарения вещества; а - радиус частицы.
Уравнение (1) состоит из пяти слагаемых, каждое из которых описывает определенный процесс, проходящий в частице под действием МЛИ. Падение на частицу излучения приводит к ее нагреву, первое слагаемое учитывает поглощение энергии в единицу времени. В результате этого вступают в силу следующие процессы: прежде всего отдача тепла в окружающую среду за счет теплопроводности, далее затраты внутренней энергии на переизлучение. Если условия позволяют частице загореться при достижении температуры воспламенения, то это учитывается в следующем
Компьютерное моделирование процесса горения.
члене. И последнее, возможно испарение частицы, но это также требует определенных условий, например, для сажистой частицы температура должна быть выше 4500 К. Уравнение (2) для радиуса также имеет несколько слагаемых. Радиус может меняться только в результате горения и испарения частицы. Прежде всего было произведено обез-размеривание системы уравнений, чтобы ЭВМ работала с числами, близкими к единице. Для этого введены новые безразмерные переменные
в = —, r = —, Qm
T, ao
QэI,
, Im = —, - = -Io t
і pRTo
?1 I о 2 ^ о 0
где д1 - теплота, выделенная в результате протекания одной из реакций горения; 10=106 Вт/м2 -характерная средняя интенсивность падающего излучения; ^ - характерное время вышорания
частицы в результате только диффузионного режима горения.
Исходная система после соответствующих преобразований примет вид
і 3
ClIm + C2Q"Ks - 0,-е2 r
-і) - C4q4 -1) - c
g3
(3)
de d-
^-= - C 6 [K s С r, в ) + K r , в ) ] (4)
d —
где Ci - безразмерные коэффициенты (i = 1, ..., 6).
Таким образом, система весьма просто выглядит в безразмерном виде, причем при этом коэффициенты, стоящие у переменных, наглядно показывают, какие процессы будут преобладать на различных стадиях. В качестве языка программирования был выбран Паскаль, метода решения дифференциальных уравнений -одна из схем Рунге-Кутта. Использовался метод Рун-ге-Кутта четвертого порядка в модификации Мерсона [3], который позволяет оценивать погрешность вычисления на каждом шаге и принимать решение об изменении его в нужную сторону, а также возможность выводить результаты вычисления с заданным шагом. Вообще схемы Рунге-Кутта очень удобны для расчетов на ЭВМ, причем все они обладают высокой точностью, являются явными и допускают расчеты с переменным шагом, что позволяет значительно повысить скорость вычисления. На примерах систем уравнений с аналитически найденными решениями тесты показали, что модификация Мерсона работает не только бъстрвй, что весьма ощутимо, но и точней при одних и тех же начальных шагах вычисления. Найден оптимальный начальный шаг, который обеспечивает высокую точность вычисления и минимальные затраты по времени. Система также была проверена на устойчивость решения, никаких отклонений замечено не было. Первоначально программа работала под операционной системой DOS, интерфейс бъш для нее написан на Turbo Vision М.Ф. Кондрашиным. Создание интерфейса под операционную систему Windows 95 на Delphi 3 позволит в значительной степени расширить возможности
данного программного продукта. Программа оформлена в виде нескольких блоков: расчетный математический и интерфейс, обеспечивающий ввод и вывод результатов. Основная часть программы - это интерфейс, который уже в свою очередь вызывает соответствующие процедуры для выполнения необходимой работы по командам пользователя. Программа в самом начале выводит характерные величины процесса, время и температуру.
Эти величины сразу же пересчитыгваются и выводятся на экран после изменения какого-либо параметра системы, что дает возможность уже по начальным условиям, не решая системы уравнений, судить о возможных температурах и временах протекания процесса. Есть еще одна особенность работы программы. Решение системы ищется на интервале от нуля до характерного времени вышорания частицы. В результате получаем зависимости температуры и радиуса на всем временном интервале протекания процесса горения, т.е. до равновесного состояния частицы и среда.
Таким образом, ничего не зная о поведении частицы, сразу получаем ответ на вопрос о реальном времени протекания процесса. При этом всегда есть возможность рассмотреть процесс на меньших или больших временных интервалах. Одним из достоинств данной программы является то, что она строит нормированный, маштабированный график зависимостей температуры и радиуса частицы от времени, где сообщается дополнительная информация о процессе и приводятся начальные условия. Дополнительно можно построить графики зависимостей от времени и температуры основных частей системы уравнений, что дает массу информации о процессах, проходящих на поверхности частицы, а именно: в каком режиме она горит, в какой момент времени она горит, а в какой нет, как зависит величина потока кислорода и вещества частицы от температуры:, как она менялась в ходе процесса и др. Если график, построенный программой, не устраивает, то можно обратиться к программам, специально предназначенным для построения графиков и обработки данных. Это возможно благодаря тому, что создается текстовым файл, в котором в табличном виде находится решение системы уравнений. В первом столбике - время, во втором - температура, в третьем - радиус. В расчетах использовались частицы с размерами от одного микрометра до одного миллиметра. При больших размерах будут возникать эффекты, не описываемые данной моделью, такие как неравномерность распределения тепла, циркуляция воздушных масс и другие. Ограничения, накладываемые на интенсивность, связаны с тем, что при плотностях мощности, больших 1010 Ватт на метр квадратный, возникает оптический
Юі
пробой, также не описываемый данной моделью.
Для расчетов была выбрана углеродная частица, так как дли нее есть наработанная экспериментальна база, позволяющая проверить теоретические расчеты. На рисунках 1, 2 представлены результаты расчета динамики температуры и радиуса частицы при различных начальных радиусах, где начальная температура частицы и температура среды Т0 = 300 К, интенсивность падающего излучения I = 2Ч108Вт/м2.
Рис. 1. Зависимость температуры частиц от времени при I = 2х108Вт/м2. Кривая 1 соответствует а(0)=100 мкм, 2 — 80 мкм, 3 — 50 мкм, 4 — 20 мкм, 5 — 8 мкм
Рис. 2. Зависимость радиуса частиц от времени при I = 2х108Вт/м2. Кривая 1 соответствует а(0)=100 мкм, 2 - 80 мкм, 3 — 50 мкм, 4 — 20 мкм, 5 — 8 мкм
Еидао, что при небольших интенсивностях изучения изменение радиуса от времени отсутствует, так как при таких условиях нет процесса горения. Это означает, что радиус частицы не мешется, кроме как в результате теплового расширения, но мы его не учитываем из-за его незначительности:. Небольшой нагрев вещества не приводил? к возгоранию, это должно проис-
ходить по достижению некоторой определенной температуры воспламенения. Диффузионный режим горен^ определяется скоростью диффузии кислорода, окислителя к поверхности облучаемой частицы, и когда она начинает тратить кислород в процессе горен^ так быстро, что он не успевает поступать к поверхности, то наступает кинетический режим горения, который характеризуется скоростью химических реакций вышорания углерода ( СО2 + С = 2СО ) . По достижению еще более высокой температуры частица тухнет, так как идет процесс активного испарения, и поток испаренного вещества не дает проникнуть окислителю к поверхности частицы.
Из рисунков видно, что прогрев частицы до максимальной температуры зависит от начального радиуса и интенсивности падающего излучения. После достижения максимальной температуры она выжодит на стационарный режим, которым характеризуется слабой зависимостью от времени. В этот период устанавливаются потоки испаренного углерода и приходящего кислорода. В данном случае процесс активного испарения идет до тех пор, пока температура поверхности частицы больше 4500 К. Далее поток углерода от частицы уменьшается и дает приток кислорода, что влечет за собой воспламенение частицы и замедление падения температуры. После чего радиус частицы уменьшается настолько, что частица расходует тепла больше, чем получает, при этом температура снова выкодит на стационарным режимі. Этой температуры уже недостаточно для интенсивного горения, а тем более испарения, поэтому далее процессы идут крайне медленно и очень долго по сравнению с интенсивным процессом горения и испарения. Основным результатом проделанной работы является программным продукт, позволяющий сравнительно просто провести расчеты основных параметров частицы в поле МЛИ и сделать соответствующие выводы о происходящих процессах в результате взаимодействия. Конкретно можно рассчитать динамику от времени: радиуса, температуры, скорости горения и испарения, плотности потока испаренного вещества, скорости изменения радиуса и температуры, коэффициента стехиометрии и теплошд^ен^ за счет горен^, а также все те же параметры в зависимости: от температуры при заданном радиусе и интенсивности.
Литература
1 Букатый В.И., СуторихинИ.А., Краснопевцев В.Н., Шайдук А.М. Воздействие лазерного излучения на твердым аэрозоль. Барнаул, 1994.
2 Шайдук А.М. Моделирование процессов взаи-
модействия мощного лазерного излучения с дисперсными системами: Дис. ... док. физ.-мат. наук. Барнаул, 1998.
31 Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. Томск, 1992.
10 2