CC BY
18УДК 537.523+535.23
В. М. Градов
ПРОГРАММНО-МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ ДЛЯ НАУЧНЫХ ИССЛЕДОВАНИЙ СИСТЕМ
С ДОМИНИРУЮЩЕЙ РОЛЬЮ РАДИАЦИОННЫХ ПРОЦЕССОВ
Рассмотрены вопросы построения вычислительных моделей, алгоритмов их реализации и соответствующего программного и информационного обеспечений для исследования процессов в системах с объемными источниками селективного излучения, использующих разрядную плазму с разными видами неравновесности. Предложен самосогласованный алгоритм расчета характеристик элементов систем, находящихся в едином радиационном поле сложного дискретно-непрерывного спектрального состава. Приведены результаты численных экспериментов, иллюстрирующие возможности вычислительного комплекса.
E-mail: gradov@bmstu.ru
Ключевые слова: программный комплекс, имитационная и аналитическая вычислительные модели, алгоритм и метод расчета, плазма, радиационный перенос, осветительная система, импульсные и дуговые лампы.
Системы, рассматриваемые в работе, могут быть идентифицированы приборами квантовой электроники, облучательными и светотехническими установками, имитаторами излучений различной физической природы, установками радиационного нагрева и др. В качестве источников селективного излучения используются газоразрядные источники лампового типа.
Физика процессов в установках на базе разрядной селективно излучающей низкотемпературной плазмы достаточно сложна, особенно если учесть, что источник излучения в составе устройства подвержен воздействию собственной многократно отраженной радиации, и его характеристики могут заметно изменяться по сравнению с работой вне системы. Влияние системы на параметры разряда настолько велико, что при оценке эффективности излучения источника оказывается затруднительным делать какие-либо прогнозы на основе одного только сопоставления спектральной кривой испускания открытого источника и спектра поглощения приемника. В этой ситуации трудно переоценить роль компьютерного эксперимента.
Построению методов расчета наиболее сложного представителя рассматриваемых систем — систем накачки лазеров — в различных приближениях посвящено большое число публикаций [1-8]. В данной работе формулируется замкнутый комплекс реалистичных нелинейных математических моделей физических процессов в системах
достаточно широкого класса, рассматриваются различные вычислительные аспекты их реализации и вопросы построения компьютерных средств расчета, ориентированных на проведение полномасштабных вычислительных экспериментов в наиболее общем случае нестационарной, неоднородной и неравновесной плазмы разрядов.
В результате появляется возможность исследования плазмы в качестве эффективного трансформатора излучения из одних спектральных компонент в другие и глубокой оптимизации систем.
Общая характеристика рассматриваемых технических устройств. Устройства представляют собой совокупность одного или нескольких разрядных источников излучения и светотехнического оборудования, обеспечивающего концентрацию лучистой энергии на приемнике излучения, например, активных элементах лазеров, нагреваемых поверхностях в оптических печах и т.д. [1-3].
Системы обладают в общем случае достаточно сложной конфигурацией и включают многочисленные элементы (рис. 1): отражатели и контротражатели, форма поверхности которых может быть самой прихотливой; различные среды (хдадагенты, иммерсию, моноблоки, оболочки источников излучения и т.д.); приемники (полупрозрачные среды, поверхности). Материалы данных элементов имеют оптические характеристики, существенно зависящие от длины волны излучения, поверхности элементов могут быть отражающими, преломляющими или поглощающими, среды — излучающими, поглощающими и рассеивающими, а отражение имеет диффузный, зеркальный или смешанный диффузно-зеркальный характер.
Источник излучения состоит из разрядной трубки, заполненной плазмообразующей многокомпонентной средой, и тракта жидкостного
Рис. 1. Система накачки лазера на АИГ:Ш:
1 — тракт охлаждения и иммерсионная среда (жидкость); 2 — оболочка лампы (кварц); 3 — плазма разряда; 4 — моноблок (кварц); 5 — отражатель (диффузный); 6 — активный элемент из алюмоиттриевого граната
или воздушного охлаждения (или без него). Оболочки, стабилизирующие разряд и образующие канал охлаждения, в результате разогрева могут стать излучающими. Материал оболочек (кварц, лейкосапфир, поликор) обладает свойствами селективного поглощения оптического излучения, а также его рассеяния (поликор). Для решения задач пассивной или активной фильтрации спектральных компонент на трубки может напыляться отражающее или поглощающее покрытие.
Центральное место в системах занимает разрядная плазма, которая служит генератором мощных лучистых потоков с КПД преобразования электрической энергии в лучистую, достигающим в ряде случаев 70—80 и даже 90%. При этом излучение распределено в широком спектральном диапазоне от ультрафиолетовой до далекой инфракрасной области и имеет сложный дискретно-непрерывный и спектральный состав. При описании спектров источников излучения необходима тщательная проработка отдельных спектральных компонент по частоте, что резко увеличивает объем вычислительных работ, так как появляется новая координата, характерная для радиационного переноса — частота излучения.
Процессы в разрядах во многих случаях являются нестационарными, плазма неоднородна по объему и в ней возможно проявление различных эффектов неравновесности. Плазма разрядов высокого и среднего давлений характеризуется давлениями 0,5 ... 4,0 МПа и температурами на оси (7... 20) • 103 К (Н^, Хе, Кг, Аг плазмообразующие среды), (0,01... 0,5) МПа и (3 ... 10) • 103 К (пары щелочных металлов Ыа, К, ЯЪ, Св). В плазме низкого давления последнее составляет от нескольких десятых до сотен единиц паскалей, а температуры электронов и атомов могут отличаться более, чем на порядок. Время прохождения тока в импульсных режимах имеет порядок десятков и сотен микросекунд.
Физические процессы в разрядах высокого и среднего давлений протекают в условиях, близких к локальному термодинамическому равновесию (ЛТР). Для разрядов низкого давления характерна сильная ионизационная неравновесность, нарушение максвелловского распределения электронной компоненты по скоростям, а распределение частиц по уровням энергии оказывается весьма далеким от больцма-новского. В установлении этого распределения важная роль принадлежит не только разнообразным столкновительным процессам, но также и излучению, в том числе и возвращаемому в разряд в результате прохождения его в системе.
Отражение указанных особенностей рассматриваемых систем в математических моделях оказывается весьма непростой задачей.
Самосогласованная математическая модель процессов в системе с пространственно неоднородными излучающими средами. При
построении комплексной модели физических процессов, протекающих в едином радиационном поле системы, надо решить прежде всего вопрос описания распространения излучения в системе. В практике решения задач радиационного теплообмена в замкнутых системах поверхностей это делают обычно на основе обобщенного зонального метода и его модификаций [4]. Существенное усложнение задачи возникает, когда, как в данном случае, в качестве источников излучения используется плазма разрядов с селективным спектром и объемным характером излучения. Тогда используют такие методы, как аналитический, графоаналитический, на основе аппарата обобщенных и локальных угловых коэффициентов, а также алгоритмы, базирующиеся на методе Монте-Карло [5]. При этом принимают ряд допущений, основанных на представлениях о нитевидных источниках и абсолютно поглощающих приемниках, приближении изотропного источника излучения с абсолютно поглощающей плазмой, зачастую ограничиваются монохроматической оценкой оптической эффективности без расчета детального энергетического баланса в системе. Кроме того, методология указанных подходов предполагает известное спектральное распределение излучения источника, которое берут по данным экспериментов для открытых излучателей, хотя в составе системы характеристики плазмы существенно изменяются. Перечисленные методы не позволяют провести полный анализ энергетики систем с частично прозрачными объемными излучателями, характеризующимися сложным спектром лучистых потоков, который вдобавок сильно изменяется в течение импульса тока.
Для комплексного расчета установок необходимо построение универсальных моделей и алгоритмов, начало которому было положено в [6-7] при разработке самосогласованного метода расчета систем накачки. При этом моделируется замкнутый круг процессов, протекающих в разряде, приемнике излучения и других элементах системы в едином поле излучения сложного спектрального состава.
Абстрагируясь от конкретного конструктивного назначения элементов, светооптическую систему сложной геометрии удобно представить в виде совокупности подобластей, содержащих ряд излучающих или поглощающих сред, которые определяются спектральными коэффициентами поглощения, преломления и усиления, а в случае плазмы — еще и термодинамическими, транспортными, теплофизически-ми и энергетическими характеристиками. Подобласти ограничиваются поверхностями первого или второго порядков, которые могут быть отражающими (зеркально, диффузно) или преломляющими. В этой модели система присутствует в виде своего оптического аналога [6].
Замкнутая модель системы строится как комбинация имитационной и аналитической моделей процессов. Каждая из излучающих подобластей, рассматриваемая как пространственно неоднородная среда,
покрывается сеткой, узлы которой являются центрами испускания. Задаются распределения параметров, определяющих излучение среды (для разряда в случае, например, справедливости приближения ЛТР — это температурное поле при фиксированных давлении и соотношении компонент). В узлах строится веер равномерно распределенных по сфере лучей, вдоль которых в соответствии с локальными характеристиками среды выпускаются группы так называемых фотонов с последующим слежением за их прохождением по системе. Формально фотон (или квант энергии) {Ф} = {АЕи,Х,Р} представляет собой выделенную порцию лучистой энергии, обладающую частотой V, энергетическим весом АЕи, вектором положения в пространстве X и направлением движения Х. Слежение за каждым фотоном продолжается до его полного поглощения в системе или ухода из нее. Фотон испускается одной из излучающих подобластей, а поглощается или усиливается всеми остальными. Фиксируя место поглощения испущенного кванта, можно найти распределение лучистой энергии по элементам системы.
Затем выполняется расчет параметров разряда и других элементов системы и вновь повторяется розыгрыш фотонов. Данная процедура в принципе может быть доведена до конца. Однако на практике такой метод расчета оказывается чрезвычайно трудоемким в свете высоких требований к точности определения дивергенции лучистого потока, которая используется затем в итерациях при расчете температурного поля. Более точно и экономично расчет переноса излучения в отдельных областях со сравнительно простой геометрией можно провести, применяя те или иные численные методы решения точного уравнения переноса излучения. При такой схеме результат моделирования прохождения совокупности фотонов в системе используется только для определения результирующего потока излучения в граничном условии для уравнения переноса в излучающих областях.
Для решения задачи лучистого переноса в плазме требуется построение моделей состояния плазмы с учетом того, что состояние плазмы зависит от радиационного поля, которое порождается самой плазмой и другими элементами системы, т.е. необходимо решать самосогласованную задачу.
Учитывая сложность задачи при практической реализации замкнутой модели системы, рассматриваются два способа согласования моделей среды и моделей радиационного переноса. Согласно первому способу, моделируются системы, для которых необходимо настраивать пространственные распределения характеристик разрядного объема. Это имеет место при исследовании систем, во-первых, в нестационарных режимах работы источников излучения; во-вторых, в случае,
когда плазма неравновесна, и, в-третьих, когда плазма находится в состоянии ЛТР и стационарна, но является оптически толстой для наиболее значимых участков спектра с точки зрения энергетики разряда и чувствительности приемника (например, разряды высокого давления в парах ртути). Сюда же относится моделирование излучения плазмы с большой ролью самообращенных резонансных линий (разряды в парах щелочных металлов).
В противоположность описанному варианту можно достичь значительного упрощения, если пренебречь пространственным распределением величин в ЛТР-плазме (разряды с относительно малой оптической плотностью в широком спектральном диапазоне, например, разряды в Хе, Кг, Аг). При этом модель строится в изотермическом приближении [6]. Задают рабочее давление и температуру плазмы, однородные по объему. После розыгрыша фотонов электрическую мощность, вкладываемую в разряд, определяют как сумму мощностей излучения, поглощенного во всех средах, кроме плазмы. Потери энергии на оболочках за счет кондуктивно-конвективного переноса теплоты устанавливают из расчетов открытого источника. По найденной плотности поглощенной мощности в элементах системы решают задачу определения температурных полей в них, а в системах накачки может быть еще и проведено моделирование работы лазера в режимах свободной генерации, запасания энергии, генерации моноимпульса и т.д.
Для реализации изложенной концепции расчета был разработан комплекс вычислительных моделей (рис. 2). Задавая параметры взаимодействия моделей, можно использовать их автономно.
Для реального вычислительного эксперимента требуется весьма полная и надежная база данных, содержащая широкий круг характеристик элементов систем и используемых материалов. Центральное место в базе занимает подсистема, связанная с информационным обеспечением моделирования процессов в плазме (рис. 3).
Опыт работы показывает, что при построении данной подсистемы необходимо обеспечивать получение всех основных характеристик плазмы хорошо апробированными методами и средствами моделирования. Связано это с тем, что, во-первых, для задач моделирования требуется точная привязка характеристик к параметрам плазмы. Имеется явный дефицит данных такого рода. Во-вторых, важным является условие полноты данных для нужного диапазона параметров. В-третьих, при поисках перспективных плазмообразующих сред такое построение базы дает мощное средство научного исследования новых проектов.
Некоторые свойства плазмы представляют собой функцию, по крайней мере, нескольких переменных. Например, коэффициент оптического поглощения для плазмы, находящейся в состоянии ЛТР,
МОДЕЛИ ПРОЦЕССОВ В СИСТЕМЕ
1 1
| ИНФОРМАЦИОННОЕ ОБЕСПЕЧЕНИЕ
БАЗА ДАННЫХ ПО КВАРЦУ И САПФИРУ
БАЗА ДАННЫХ ПО СВОЙСТВАМ ЭЛЕМЕНТОВ
БАЗА ДАННЫХ ПО СВОЙСТВАМ
АКТИВНЫХ СРЕД _
Рис. 2. Структура разработанных математических моделей физических процессов в рассматриваемых системах
Рис. 3. Структура разделов программно-математического обеспечения информационной поддержки моделирования плазмы
зависит от четырех переменных, причем одной из переменных является частота излучения (длина волны), по которой требуется выполнять, как правило, весьма детальную проработку спектра. Известно немного публикаций подробных таблиц данного коэффициента, например работа [9]. На рис.4 в качестве примера представлен коэффициент
1.00Е+06 г
^ 1.00Е+05 -оГ
| 1.00Е+04 -си
3"
§ 1.00Е+03 -о
н 1.00Е+02 -0)
| 1.00Е+01 --8-
1.00Е+00 -
о
1.00Е-01 I-0
Рис.4. Коэффициент поглощения плазмы С§-Щ-Хе. Массовые соотношения компонент плазмы и С§:Хе 1:1 и 1:1,5 соответственно. Давление плазмы
р = 0,5 МПа, температура Т = 4000 К
1 1
1 Ii
1 Ду ЛА
у v w V УЧ
400 800 1200 1600 2000
Длина волны, нм
поглощения с детализацией около 1500 узлов по частоте, с которой он хранится в базе и которая позволяет адекватно учесть дискретную составляющую спектра. Наконец, для неравновесной плазмы большинство свойств плазменной среды необходимо рассчитывать в самом процессе моделирования (например, тот же коэффициент поглощения, коэффициенты диффузии, скорости неупругих процессов и т. д.), и эти свойства становятся неотъемлемой частью модели (см. рис. 2).
Помимо подсистемы, обслуживающей моделирование разрядов, необходимы также другие подсистемы, поставляющие информацию о радиационных и теплофизических свойствах материалов и поверхностей элементов устройств. Сюда относятся спектральные коэффициенты поглощения, преломления и рассеяния материалов при различных температурах, спектральные коэффициенты поглощения, преломления, отражения и степени черноты поверхностей с учетом их состояния и т.д. Эту часть общей базы данных допустимо пополнять и обновлять из различных литературных источников.
Перенос излучения сложного спектрального состава в излучающих средах произвольной оптической плотности. Лучистый перенос в плазме разрядов и других излучающих подобластях рассматриваемых систем описывается на основе уравнения переноса. В общем случае поглощающих, излучающих и рассеивающих сред данное уравнение в приближении некогерентного рассеяния для частоты излучения V на основании [4, 10] может быть записано в следующем виде:
^^^^ + (Х) + ^ (Х)) IV (Х, 8) =
с2 1 Г
= ^(Х)(1 + ^IV(Х, 8)) + —^(Х)J Рр(Х, 8, 8%Г, (1)
4п
где к1,(Х),^(Г) — спектральные коэффициенты поглощения и рассеяния среды; г — радиус-вектор точки; Х — единичный вектор направления распространения излучения; рр(Х, Х, Х') — индикатриса рассеяния с направления Х на направление Х; ¿П' — элементарный телесный угол; jV (Г) — объемный коэффициент собственного излучения среды; с — скорость света. Сомножитель при jV (Г) отвечает за индуцированное (вынужденное) испускание.
При построении моделей учитывается, что в селективно излучающей плазме значительная роль в переносе излучения принадлежит спектральным линиям, которые уширены различными механизмами (штарковским - электронами и ионами, резонансным, вандерваальсов-ским, доплеровским). Перенос линейчатой составляющей излучения
происходит на непрерывном фоне, обусловленном эффектами фоторекомбинации электронов и ионов, а также тормозными процессами в полях нейтральных частиц и ионов.
Из уравнения (1) в результате интегрирования по всем направлениям (телесному углу) получают уравнение сохранения энергии излучения
где UV = ^^П — объемная плотность энергии излучения; FV —
плотность спектрального потока излучения; к'1/(Г) — суммарный коэффициент поглощения плазмы на частоте V, исправленный на вынужденное излучение.
Учет переноса излучения в случае реальных спектров даже при наличии ЛТР приводит к ряду вычислительных трудностей, возрастающих при расчетах неплоских конфигураций (например цилиндрических). Эти трудности связаны с тем, что точное уравнение переноса излучения записывается для спектральной интенсивности излучения на фиксированной частоте в выделенном направлении в пространстве. В то же время в уравнении энергии для электронов плазмы фигурирует дивергенция интегрального по спектру лучистого потока, поэтому требуется проводить интегрирование уравнения переноса и по пространству (углам), и по частоте. В настоящей работе проблема интегрирования по пространству решается путем непосредственного взятия многомерных интегралов по углам в каждой точке излучающе-поглощающей среды или сведением задачи к дифференциальным уравнениям — диффузионному приближению или приближению Шустера-Шварцшильда [4, 11, 12]. Вопрос интегрирования по частоте возникает с дискретным спектром и решается путем индивидуального учета наиболее интенсивных линий, а большое число слабых линий "размазывается" по занимаемым ими частотным интервалам [5, 13, 14].
Интенсивность излучения в каждой радиальной точке разрядного объема зависит от сферических угловых координат в и ее определяют по формуле
divF, = 4njv (r) - ckV (r)Uv,
(2)
l(r,e,ip) / l(R,9,ip)
+
где 1-(Я, в, ф) — интенсивность внешнего излучения на границе плазменного объема, поступающего внутрь в направлении, которое определяется углами вф, (для открытого источника излучения эта величина равна 0); 1(т, в, ф) — длина луча от точки с радиальной координатой т до поверхности цилиндра в направлении вф; в, в'- координаты вдоль луча.
Система уравнений лучистого переноса в диффузионном приближении включает уравнение (2) и уравнение для плотности потока, которое в общем случае записывается в виде [10]
Ру = — ^, (4)
V 3к^ Ат '
п
где к^ = к5и + 8иви; к^и 8и = 1 - 0,27(в)вт(2в)Ав — ко-
0
эффициенты оптического поглощения с учетом рассеяния, скорость света, коэффициенты рассеяния и формы индикатрисы рассеяния [4], соответственно.
Граничные условия к модели (2), (4) при расчетах разрядов вне
отражающих систем имеют вид:
0 0
т = 0, —¡— = 0,
ат
„ тт А 1 + Ри ¿и
т = л,ии = —— •
к 1 — Ат'
где pV - спектральный коэффициент диффузного отражения на границе плазма-оболочка; А — константа, оптимизация которой позволяет повысить точность диффузионного приближения.
При моделировании процессов в плазме в составе системы, когда имеет место возврат лучистой энергии в разряд, граничные условия при т = Я формулируются иначе:
с аи
т = Я, • -и = — ^р = ^с(1 — фи), (5)
3к' ат
где Рис, — поверхностные плотности спектральных потоков собственного излучения плазмы и излучения, поглощенного в разряде в
I Vр
процессе многократных проходов излучения в системе; = —--
р vc
коэффициент возврата излучения в плазму разряда.
При моделировании кинетики заселения уровней энергии возникает задача расчета пленения излучения в линиях с большим коэффициентом поглощения — задача излучательного или радиационного переноса возбуждения. Эти линии имеют, как правило, в качестве нижнего уровня перехода — основной уровень (обычно резонансные
линии). Роль эффектов пленения в заселении верхних уровней таких линий чрезвычайно велика. Обычно в расчетах используют достаточно приближенные выражения для определения вероятности вылета фотона за пределы плазменного объема, основанные на решении уравнения Бибермана-Холстейна для модели однородной плазмы [15]. Такой подход для наших задач в ряде случаев приводит к заметным погрешностям, и задачу переноса радиационного возбуждения частиц приходится ставить следующим образом. Записывается стационарное уравнение неразрывности для частиц, находящихся на верхнем уровне 1 линии, образованной переходом в основное состояние атома 0, в виде
1 d , ^ dni(r). v-^ , w / ч , ч
(rDi ) + Е n*(r)(w»i(r) + Ai) —
- ni(r)(E wij(r) + Aio) + hV^y kV(r)U(r)d- = 0, (6)
j Av
где n¿(r)- концентрация частиц на уровне ¿в точке с радиальной координатой r; Di- коэффициент диффузии атомов; (r) — вероятность перехода между уровнями ¿ик за счет столкновений с частицами; kV (r) — коэффициент поглощения на частоте - в пределах рассматриваемой линии в точке r; Av — спектральная область линии.
Входящая в уравнение (6) объемная плотность энергии излучения вычисляется согласно выражению
п п/2 l(r,e,v)
UV (r) = h—0—"У d^y d^ sin ^ J ni(s)pi0(v, s)x o o o
i(r,e,ip)
x exp ^ — J kV(s')ds'^ ds. (7)
Уравнение (6) является интегро-дифференциальным уравнением, которое решается численными методами. Заметим, что во многих задачах можно пренебречь вынужденным излучением, и тогда коэффициент поглощения не зависит от населенности верхнего уровня линии Пь а записанное уравнение оказывается линейным. Различные упрощения в решении уравнения (6), в том числе и при использовании диффузионного приближения, описаны в [5].
Модель распространения монохроматического излучения в многоэлементных системах сложной геометрии. Разрядная плазма с изложенных позиций представляет собой подобласть поглощающе-излучающей среды, ограниченную цилиндрической поверхностью.
Окружность столба разряда разбивают на Иу равномерно размещенных по окружности 2пЯд точек, где д — часть разряда в азимутальном направлении, попавшая в оптический аналог.
Из каждой точки как из центра строится веер N равномерно распределенных по полусфере лучей. По каждому направлению веера выпускается Ир усредняющих квантов мощности ("фотонов") с энергетическим весом ДЕи0(й) = 1и ^П^ /Ир, где единичный вектор П задает направление вылета кванта, а интенсивность излучения на частоте V, испускаемого с площадки Д£ в телесном угле ДП в направлении П, рассчитывается согласно (3)
Ь(П) Ь(П)
1и(П) = { J 3у(х)ехр(- [ к(х')д,х')д,х\ ДйДБ,
где Ь(й) — длина луча в плазме, на котором генерируется излучение;
— коэффициент излучения плазмы в точке с координатой х вдоль луча,
Дй Дй = р (у,й )Д/, Д1 = а - а2
2N - 1 х ' 7 ' Na
2RP(Y,Ü) g _ 2пд
1 - P (Z,Ü )2
V > , —g Nv
здесь Р(У,й),Р) — проекция единичного вектора направления луча на ось симметрии полусферы У и продольную ось системы; а1,а2- координаты начала и конца системы вдоль продольной оси; Иа — число участков разбиения лампы по длине.
Для плазмы, находящейся в состоянии локального термодинамического равновесия коэффициент излучения просто выражается по закону Кирхгофа через функцию распределения Планка
З (х) = ки (х)1„р(Т (х)),
где Т(х) — распределение температуры вдоль луча; 1ир — интенсивность равновесного излучения.
После старта производится моделирование траектории фотона в системе, в ходе которого производится последовательный расчет указанного выше набора. Радиус-вектор места встречи с поверхностью вычисляется в виде
х4 = х + Р • г,
где г — расстояние от точки, определяемой вектором X до ближайшей поверхности в направлении Р.
При встрече с поверхностью фотон может отразиться, поглотиться, преломиться. В случае отражения от зеркальной поверхности новое
направление фотона вычисляется по формуле
Р! = Р — 2п(Рй),
где п — вектор нормали в точке падения фотона. При отражении от диффузной поверхности направление отражения случайно и равномерно распределено по полусфере, что алгоритмически реализуется заданием случайных чисел. В результате отражения текущий вес фотона изменяется
ДЕ^ = ,
где ДЕ^ — вес фотона до отражения; р^ — коэффициент отражения излучения от ]-й поверхности.
Пропускающая поверхность может обладать зеркальными, диффузными или диффузно-зеркальными свойствами отражения и преломления. Судьба фотона, попавшего на такую поверхность, определяется также с помощью операции со случайными числами: если очередное случайное число 7, генерируемое датчиком случайных чисел, удовлетворяет условию 7 ^ рVj, то фотон отражается, иначе — преломляется в другую среду. В случае зеркальной поверхности коэффициент отражения рассчитывается по формулам Френеля, а иначе берется из литературных источников. Если граница раздела диффузно-пропускающая, то направления отраженного и преломленного фотонов случайны и равномерно распределены по полусферам отражения или преломления.
Текущий вес фотона после прохода расстояния £ в ]-й среде с коэффициентом поглощения ^ составляет
ДEVj = ДEVj-! ехр — / ^ (X
I. 'о
Слежение за фотоном продолжается до тех пор, пока его текущий вес не уменьшится до задаваемого пользователем значения. После того, как все выпущенные фотоны распределятся по элементам системы, будет найдена интегральная по спектру дивергенция лучистого потока, тогда проводится расчет температурного поля в среде.
Модель неравновесной пространственно неоднородной многокомпонентной плазмы. В настоящей работе модель неравновесной излучающей плазмы строится в гидродинамическом приближении [16], т.е. плазма представляется в виде совокупности газов, которая в типичном для исследуемых устройств варианте включает следующие компоненты: электроны, ионы, атомы (в основном состоянии и возбужденные на различные уровни энергетического спектра) и нейтральный газ-буфер. Данное описание оправдано, если частота столкновений частиц достаточно велика. Компоненты находятся во
взаимодействии друг с другом, испытывая упругие и неупругие столкновения. Через плазму идет мощный поток энергии с трансформацией в различные ее виды.
Уравнение непрерывности для частиц одного сорта а записывают в виде
дп -*
+ &у[па(м + иа)] = Па,
где па — скорость образования (или распада) частиц сорта а в единице объема в результате неупругих столкновений с другими сортами частиц; и — среднемассовая скорость движения плазмы; иа — диффузионная скорость частиц.
Уравнения переноса импульса без учета вязкости и термодиффузии имеют вид:
Ши
Ш
ШеПе^- = -en-gE - Vpe - mene Е Vek (Ue - Uk); (8)
k=e
ш^ — = en^ - Vpi - n E ШгкVik(Ui - Uk); (9)
du "dt"
k=i
Шп nn ^Uf = -Vpn - n^E mnk Vnk (Un - Uk ), (10)
k=n
здесь , mjk, nj, е, E,pj, Р^ — соответственно, масса частицы сорта приведенная масса частиц сортов ] и к, объемная концентрация частиц -го вида, заряд электрона, напряженность электрического поля, парциальное давление ] -й компоненты, средняя частота упругих столкновений частиц ] и к с передачей импульса; индексы е, г, п — электроны, ионы и нейтральные частицы, возбужденные до п-го уровня энергии. Суммирование в правой части выполняется по всем сортам частиц, с которыми взаимодействуют частицы рассматриваемого сорта. Суб-
Ш д ^
станциональная (лагранжева) производная — = —- + и V.
ш: д£
Складывая уравнения (8) и (9) и опуская при этом инерционные члены, можно получить выражение для радиального потока амбипо-лярной диффузии заряженных частиц к непроводящим стенкам оболочки, стабилизирующей разряд. Полные уравнения сохранения для плазмы как целого получают суммированием уравнений (8)-(10) по всем компонентам.
При записи уравнений энергии для легких и тяжелых частиц под температурой будем понимать величину, которая после умножения на константу Больцмана дает две трети кинетической энергии хаотического движения частиц. Для электронного газа уравнение энергии с учетом неупругих потерь, связанных с возбуждением тяжелых частиц
и излучением, можно представить в виде [16]
d dt
ne
5 kTe+I^j divU = -div [Fe+{5 kTe+I)neUe
2 kTe +l)] +n
oE2 - ve^T 3 k(Te - Th)2—Vek -° mk
+
^ 2
k
Emn m divFR. (11)
Уравнение энергии для тяжелых частиц
d ( ^ \ 5IrTl1. ^ d^ynh 2 kThj + nh 2 kThdivu =
= -divFh + ne V 3к(Te - Th)2—vek, (12)
и 2 —k
k
здесь nh = n¡ — суммарная концентрация всех тяжелых частиц.
i
Суммирование в правой части выполняется по всем сортам тяжелых частиц.
Для вычисления скоростей образования или гибели частиц naв правой части уравнения непрерывности необходимо строить модель кинетики заселения возбужденных уровней энергии излучающей компоненты [15, 16]. Наиболее разработана модель ударно-излучательной ионизации и рекомбинации в модифицированном диффузионном приближении [15]. В задачах моделирования излучения в отдельных спектральных линиях, которые решаются в настоящей работе, модель кинетики приходится строить на основе индивидуального учета уровней, особенно низкорасположенных. При этом большое число уровней, примыкающих к континууму, объединяются в блок со свободными электронами.
При расчете скоростей процессов возбуждения и ионизации в результате столкновений со свободными электронами необходимо вычислять интеграл свертки сечений неупругих процессов и функции распределения электронов по скоростям. Высокоэнергетичный "хвост" этой функции может заметно отличаться от равновесного. На рис. 5 видно, что в области энергий возбуждения резонансного уровня (примерно 5 эВ) функция распределения резко отклоняется от максвел-ловской, при энергии электронов 10 эВ различие доходит до 70. В настоящей работе изотропную часть функции распределения электронов находят согласно [16] из уравнения
de
df 0
ea(e) f (e) + b(e)f 0(e)
= c(e)
f 0(e) - ^ f 0(e - eoi) glVo
(13)
0.8
Ä6 J
Ш 0.2
1 X % \ X
----2 о
* \ * \ \ \ \
\ \ \ \ \ \
\ \ \ * ч ч ч
л \ > * ч
0 5 10 15 20 25 „ 30
е, эВ
Рис. 5. Отношение истинной функции распределения электронов по энергиям к максвелловской на оси разрядной трубки радиусом Я = 9,5 мм в разряде в парах Щ с Аг. Давление Щ и Аг, соответственно, 0,85 и 465 Па; значения тока
1: 1 — 0,6 А; 2 - 0,4 А; 3 — 0,2 А
где а(г), Ь(г) — коэффициенты, которые являются функционалами функции распределения; г — энергия электронов; п0 , п- концентрации частиц излучающей компоненты на основном и резонансном уровнях; д0,0! — статистические веса основного и резонансного уровней, с(г) — см. [16].
На основе указанных групп уравнений движения, энергии, непрерывности, кинетики возбуждения и ионизации, переноса излучения строятся математические модели процессов в разрядах с различными характеристиками неравновесности плазмы. В общем случае все фигурирующие в модели материальные функции приходится рассчитывать в ходе решения задачи. Предварительное представление величин в виде таблиц практически не может быть выполнена.
Рассмотренная математическая модель неравновесного разряда носит достаточно общий характер и может быть использована при проведении вычислительных экспериментов для исследования источников излучения различных типов в широком диапазоне рабочих параметров.
На основе данной модели могут быть построены более простые модели, в которых ряд процессов описывается с применением менее сложных схем. Так, большой класс разрядных источников излучения высокого давления допускает описание на базе модели, в которой распределение частиц по уровням энергии вплоть до резонансного может быть принято больцмановским с температурой заселения, равной электронной. К таким источникам относятся лампы на основе Хе, Кг, ^ с рабочим давлением порядка 0,1... 5 МПа. При этом плазмообразу-ющая среда является однокомпонентной, т.е. буферный газ в разряде
отсутствует. В качестве дополнительного соотношения используется либо задание значения рабочего давления р, либо начального давления наполнения р0.
Дальнейшее упрощение модели связано с принятием больцманов-ского распределения частиц с температурой электронов по всем уровням, вплоть до основного. Неравновесность плазмы в этой модели проявляется в различии температур легких и тяжелых частиц и отклонении концентрации частиц от ионизационного равновесия вследствие диффузионных процессов. В этом случае можно предварительно составить таблицу значений коэффициента оптического поглощения, отнеся последний на одну излучающую частицу, в функции электронной температуры и концентрации электронов, что сильно экономит машинное время.
Компьютерная реализация математических моделей. Основным способом решения систем уравнений моделей являются численные методы. При этом дифференциальные и интегро-дифференциаль-ные уравнения решаются методом конечных разностей [17].
Для параболических уравнений строятся двухслойные консервативные однородные разностные схемы неявного типа. Обыкновенные дифференциальные уравнения второго порядка аппроксимируются на трехточечном шаблоне и сводятся также к системе алгебраических уравнений с трехдиагональной матрицей. Системы разностных уравнений решаются методом прогонки в различных вариантах (матричная, встречная, циклическая). Для квазилинейных разностных схем применяются итерационные процедуры разного типа (с релаксацией, фиксированным числом итераций). В ряде случаев проводится линеаризация методом Ньютона. Вопросы устойчивости и сходимости разностных решений анализируются в соответствии с [17].
Алгоритм решения систем дифференциальных уравнений строится на основе метода последовательных прогонок. В ряде случаев из разностных аппроксимаций каждого уравнения исходной системы составляется общая система алгебраических уравнений, которая решается прямым методом Гаусса, методом Гаусса с итерациями, методом Зейделя и его модификациями и др. Двумерные и трехмерные дифференциальные уравнения решаются методом переменных направлений.
Задачи о переносе излучения в приближении Шустера-Шварц-шильда и развитии процессов в разрядном контуре ставят как задачи Коши. Здесь оказываются удобными методы типа предиктор-корректор, Рунге-Кутта четвертого порядка точности, а также неявные схемы типа схем Гира. Схема Гира особенно актуальна при решении системы нелинейных электротехнических уравнений внешней цепи, так как на начальной стадии разряда из-за высокого омического сопротивления плазменного столба данная система оказывается жесткой.
При моделировании нестационарных процессов в плазме возникает вопрос о начальных условиях задачи. В случае импульсно-периодических разрядов или режимов с дежурной дугой решение данного вопроса очевидно. При одиночных импульсах начальное условие должно быть задано в общем случае исходя из состояния плазмы после пробоя разрядного промежутка и формирования проводящего канала. Эта стадия разряда в работе не рассматривается вследствие ее короткой длительности и слабого влияния на энергетику систем при длительностях более 100 мкс.
Решение отдельных уравнений математических моделей при задании в качестве параметров величин, определяемых из других уравнений, приводит после некоторых усилий, как правило, к положительному результату. Иначе обстоит дело с решением систем уравнений — численной реализацией модели. В этом случае успех не гарантирован и во многом зависит от рациональной организации итерационного процесса, который представляет собой множество вложенных итерационных циклов. Рассмотрим общие схемы реализации представленных выше моделей.
В модели нестационарного ЛТР-разряда в составе внешнего контура самый внешний цикл итераций осуществляется по параметрам электрической цепи. При этом настраивается омическое сопротивление Яр, определяемое в случае ЛТР-плазмы температурным профилем и давлением в плазме, и находят значения электрического тока и напряжения. Во внутренних итерациях по определению температурного поля в разряде проводят моделирование переноса излучения в системе для каждого спектрального интервала, и формируется граничное условие (5) для уравнений переноса излучения в плазме (в диффузионном приближении). Затем решаются уравнения (2), (4) переноса излучения в разряде для каждого спектрального интервала и определяется профиль дивергенции интегрального лучистого потока путем интегрирования по частоте (2), который используется при расчете профиля температур.
На рис. 6-9 представлены результаты моделирования системы накачки с одной импульсной ксеноновой лампой и двумя активными элементами АИГ: № (см. рис. 1). Радиус разрядного канала лампы составляет 7 мм, толщина стенки оболочки - 1 мм, начальное давление наполнения — 0,04 мПа, диаметр активного стержня — 6,3 мм, длина разрядного промежутка и стержней — 120 мм. В качестве охлаждающей жидкости и иммерсионной среды рассматривается дистиллированная вода. Отражатель принят диффузным, его спектральный коэффициент отражения соответствует засыпке из диоксида циркония. Разрядный контур имеет следующие параметры: емкость конденсатора С = 200 мкФ, индуктивность контура Ь = 20 мкГн, начальное напряжение на конденсаторе и0 = 1200 В. Лампа работает с дежурной
1.40Е+04
1.20Е+04
1.00Е+04
8.00Е+03
6.00Е+03
4.00Е+03
0.0 0.2 0.4 0.6 0.8 1.0
г
Рис. 6. Температурное поле в разряде в максимуме разрядного тока (момент времени £ = 90 мкс). Кривые 1-4 соответствуют номерам вариантов параметров системы накачки, приведенным в таблице
дугой в режиме периодического следования импульсов. Спектральные коэффициенты поглощения кварца типа КУ-1, активного и неактивного поглощения в активной среде и воде, отражения отражателя взяты из работ [18, 19]. В расчетах варьировались оптические свойства материалов элементов системы, при этом соответственно изменялись коэффициенты фи возврата излучения в плазму разряда. Усредненные по областям спектра и времени импульса значения для четырех различных вариантов выполненных расчетов приведены в таблице.
Таблица
Номер варианта Длина волны, нм
100-360 360-480 480-620 620-712 712-917 917-3000 3000-6000
1 0 0 0 0 0 0 0
2 0,043 0,43 0,09 0,42 0,10 0,40 0
3 0,04 0,72 0,14 0,70 0,16 0,69 0
4 0,038 0,950 0,18 0,91 0,20 0,90 0
Значение = 0 соответствует варианту открытой лампы (вне
......
* V \\\
-1 -----2 3 ----4 \ V г и
\ \
да
Рис. 7. Спектральное распределение энергии излучения лампы в составе системы накачки. Кривые 1-4 соответствуют номерам вариантов параметров системы накачки, приведенным в таблице
системы). В основном возврат излучения происходит в областях вне полос поглощения приемника (активной среды). В коротковолновой области кварц и активная среда имеют высокий коэффициент поглощения и здесь « 0. То же имеет место для кварца в инфракрасной области. Чем прозрачнее материалы системы накачки и выше коэффициент отражения отражателя, тем больше значение , т.е. тем лучше условия возврата излучения в разряд. Отметим, что с точки зрения эффективности системы накачки высокие значения коэффициента в области полос поглощения АИГ:№ свидетельствуют о неблагоприятной ситуации (низкое качество системы накачки). В то же время данный факт положителен для остальных (непроизводительных) областей спектра.
На рис. 6 видно, что возврат излучения в разряд заметно изменяет профиль и уровень температур. В результате происходит трансформация спектра в коротковолновую сторону (см. рис. 7): в области около 200... 400 нм энергия, теряемая разрядом, возрастает почти в 2 раза, несколько меньше увеличивается излучение в полосах поглощения неодима. Следует отметить, что при этом повышаются энергетические потери на оболочке, т.е. соответствует известному из экспериментов факту снижения предельной энергии ламп в составе систем. На рис. 8
15000
Рис. 8. Зависимость температурного профиля в разряде от радиальной координаты в течение импульса тока. Данные соответствуют варианту 2, приведенному в таблице
t, с Л.,пгп
Рис. 9. Изменение спектрального распределения результирующего потока излучения на внутренней поверхности разрядной трубки в течение первой половины импульса тока (до 100 мкс) в зависимости от времени. Данные соответствуют варианту 2, приведенному в таблице
и 9 показаны временные развертки температурного поля и спектрального распределения излучения. Данные этих рисунков подтверждают возможности вычислительного эксперимента при исследовании ряда характеристик плазмы, которые очень сложно измерить экспериментально.
В ситуации, когда можно проводить расчеты в изотермическом приближении (разряды в инертных газах), вся вычислительная процедура заметно упрощается, так как уравнение энергии превращается в обыкновенное дифференциальное уравнение и включается в общую схему решения задачи Коши; параметры, связанные с распространением излучения в системе, рассчитывают заранее при фиксированных температурах и давлении в плазме и их используют в модели в виде массивов данных.
При решении задач, в которых рассматривают процессы в неравновесной плазме, вычисления заметно усложняются. Применительно к стационарным открытым разрядам (вне системы) внешний цикл итераций организуется по концентрации электронов на стенке, которую для начала процесса задают в известной мере произвольно. Такое построение алгоритма продиктовано сложностью численного решения уравнений, описывающих амбиполярную диффузию заряженных частиц при краевых условиях второго и третьего рода. В ходе расчета концентрация электронов на стенке настраивается в соответствии с граничным условием задачи. Следующий вложенный цикл итераций проводится по температуре электронов. Дальнейшее вложение циклов производится по следующей схеме: расчет скоростей ударных и радиационных процессов и определение поля концентрации тяжелых частиц, расчет температуры тяжелых частиц, вычисление профиля концентрации электронов, соответствующего известным к данному моменту распределениям температур легких и тяжелых частиц и концентраций тяжелых частиц, определение дивергенции интегрального по спектру лучистого потока из уравнений переноса излучения. Выход из итераций по температурам компонент и концентрации электронов осуществляется по балансу потоков энергии и частиц в каждом узле разностной сетки.
Вся совокупность моделей (см. рис. 1) и методов расчета характеристик плазмы реализована в программно-математическом обеспечении. Организация программного комплекса выполнена в соответствии с задачами, которые приходится решать в ходе разработки и исследования различных устройств. В число этих задач входят задачи поискового характера, целью которых в основном является подбор новых плазмообразующих сред, отработка новых конструкций и новых материалов; задачи исследования процессов в существующих приборах, в том числе в критических и запредельных режимах, и, наконец, задачи
проектные, связанные с проведением серийных расчетов по оптимизации параметров и режимов работы установок.
При решении задач поискового характера требуется проводить полный цикл моделирования, начиная с вычисления параметров элементарных процессов и заканчивая определением основных выходных характеристик объекта исследования (разряда, источника излучения или прибора в целом). Задачи второй группы выполняются на существующей базе данных или путем ее расширения с учетом новых диапазонов параметров плазмы. Серийные расчеты проводятся, по возможности, по наиболее простым для каждой задачи моделям, так как требуются большие затраты машинного времени в силу своей массовости.
Структура программного комплекса такова, что в его состав входят три библиотеки программ, база данных и архив результатов моделирования.
Первая библиотека включает программы, которые могут функционировать автономно и полностью готовы к выполнению соответствующей проектной операции. В данной библиотеке представлены программы для моделирования нестационарного ЛТР-разряда в составе разрядного контура, разрядов среднего давления с частично неравновесной плазмой (в том числе с системой излучающих оболочек), разряда низкого давления с неравновесной плазмой, переноса излучения в системе с источниками излучения с неоднородным и однородным плазменным столбом, процессов в активной среде лазеров.
Вторая библиотека содержит программы для расчета термодинамических, теплофизических и оптических свойств плазмы (статистических сумм, составов в приближении ЛТР или двухтемпературном приближении, теплоемкости, плотности, коэффициентов тепло- и электропроводности, штарковских ширин линий, вероятностей спонтанного распада уровней, коэффициентов оптического поглощения в линиях и составляющих непрерывного спектра). В этой же библиотеке помещены процедуры для вычисления сечений неупругих столкновений электронов с атомами, скоростей ударных процессов, сечений фотоионизации, контуров спектральных линий, функции распределения электронов по энергиям и ряд других сервисных процедур, например, предназначенных для разбиения спектра на интервалы.
Третья библиотека содержит набор процедур для выполнения стандартных математических операций и ряда специальных вычислений (многомерная интерполяция, решение систем линейных алгебраических уравнений, одномерное и многомерное интегрирование таблично заданных функций и др.).
Результаты моделирования сохраняются в архиве. Помимо своего чисто информационного значения для специалиста в предметной
области данные из архива в ряде случаев дают возможность сделать удачный выбор начального приближения при моделировании сложных процессов, например, в неравновесных разрядах.
Понятно, что все численные методы, алгоритмы и программы тщательно проверялись путем сопоставления численных и аналитических решений, применения разных методов к исследованию одной задачи, сравнения с обширным экспериментальным материалом по всем функциональным характеристикам устройств и свойствам плазмы [5].
Вывод. Разработанный многофункциональный вычислительный комплекс служит удобным и гибким средством автоматизации научных исследований и проведения вычислительных экспериментов в области квантовой электроники, свето- и облучательной техники. Его профессиональный уровень прошел всестороннюю апробацию при выполнении большого числа разнообразных проектов в указанных ранее областях [5-8, 20].
СПИСОК ЛИТЕРАТУРЫ
1. Справочная книга по светотехнике / под ред. Ю.Б. Айзенберга. - М.: Энергоатомиздат, 1995. - 526 с.
2. С а р ы ч е в Г. С. Облучательные светотехнические установки. - М.: Энерго-атомиздат, 1992. - 108 с.
3. Лазеры на неодимовом стекле / А.А. Мак, Л.Н. Сомс, В.А. Фромзель и др. - М.: Наука, 1990.-288 с.
4. О ц и с и к М. Н. Сложный теплообмен. - М.: Мир, 1976. - 616 с.
5. Г р а д о вВ. М. Разработка методов расчета и исследование радиационных процессов в системах с разрядными источниками селективного излучения: Дис. ... д-ра техн. наук / МГТУ им. Н.Э. Баумана. - М., 2002. - 323 с.
6. Методы расчета систем накачки твердотельных лазеров / В.М. Градов, В.Г. Дорогов, Б.А. Константинов и др. // Обзоры по электронной технике: Электровакуумные и газоразрядные приборы. - 1978. - № 10 (539). - 68 с.
7. Метод самосогласованного расчета систем накачки лазеров: Характеристики систем накачки с неоднородным распределением температуры в разряде ламп / В.М. Градов, А.А. Мак, А.А. Щербаков и др. // Оптика и спектроскопия. - 1984. -Т. 57, № 1. - С. 105-109.
8. Машинное моделирование процессов в твердотельных лазерных излучателях и усилителях с ламповой накачкой: Системы накачки лазеров / В.М. Градов, Г.И. Кромский, А.А. Мак и др. - М., 1984. - 72 с. - (Препринт ФИАН, № 105).
9. Оптические свойства горячего воздуха / И.В. Авилова, Л.М. Биберман, В.С. Воробьев и др.; под ред. Л.М. Бибермана. - М.: Наука, 1970. - 320 с.
10. Зельдович Я. Б., Райзер Ю. П. Физика ударных волн и высокотемпературных гидродинамических явлений. - М.: ФИЗМАТЛИТ, 2008. - 656 с.
11. Радиационный теплоперенос в высокотемпературных газах: Справочник / И.Ф. Головнев, В.П. Замураев, С.С. Кацнельсон и др.; Под ред. Р.И. Солоухина. - М.: Энергоатомиздат, 1984. - 256 с.
12. Р е к и н А. Д. Уравнения переноса излучения в приближении Шустера - Швар-цшильда для задач со сферической и цилиндрической симметрией // ТВТ. - 1978. -Т. 16, №4.-С. 811-818.
13. Энциклопедия низкотемпературной плазмы / Под ред. В.Е. Фортова. -М.: Наука, 2001. - (Вводный том: Кн. 2). - С. 126-312.
14. С у р ж и к о в С. Т. Тепловое излучение газов и плазмы. - М: Изд-во МГТУ им. Н.Э. Баумана. - 2004. - 544 с.
15. БиберманЛ. М., Воробьев В. С., Якубов И. Т. Кинетика неравновесной низкотемпературной плазмы. - М.: Наука, 1982. - 375 с.
16. Митчнер М., К р у г е р Ч. Частично ионизованные газы: Пер. с англ. / Под ред. А.А. Иванова. - М.: Мир, 1976. - 496 с.
17. С а м а р с к и й А. А. Теория разностных схем. - М.: Наука, 2009. - 656 с.
18. Излучательные свойства твердых материалов: Справочник / под общ. ред. А.Е. Шейндлина. - М.: Энергия, 1974. - 472 с.
19. Термооптика твердотельных лазеров / А.В. Мезенов, Л.Н. Сомс, А.И. Степанов. - Л.: Машиностроение, 1986. - 199 с.
20. Математическое моделирование и исследование импульсных разрядных ламп инфракрасного излучения / С.В. Гавриш, В.М. Градов, А.В. Кузнецова и др. // Светотехника. - 2008. - № 5. - С. 14-18.
Статья поступила в редакцию 10.05.2012
