Radial luminescence structure simulation on length of the unipolar gas breakdown discharge Текст научной статьи по специальности «Физика»

УДК 537.525

Сухов Андрей Константинович

кандидат физико-математических наук, доцент Костромской государственный университет им. Н.А. Некрасова

priklmath@ksu. edu. ru

МОДЕЛИРОВАНИЕ РАДИАЛЬНОГО ПРОФИЛЯ СВЕЧЕНИЯ ПО ДЛИНЕ РАЗРЯДА УНИПОЛЯРНОГО ПРОБОЯ ГАЗА

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

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

Введение

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

Данный тип разряда, являясь безэлектродным, представляет практический интерес для создания «чистой» плазмы ввиду отсутствия металлических электродов, контактирующих с разрядом. В отличие от высокочастотных, эти разряды имеют более низкую частоту следования импульсов. Кроме того, они по своей природе очень чувствительны к внешним параметрам, что позволяет их использовать для контроля технологических плазменных процессов.

Разряд униполярного пробоя газа (УПГ) возникает в трубке с разреженным газом под действием высоковольтных импульсов потенциала одной полярности, прикладываемых к единственному внешнему покрытие-электроду (ПЭ) [1-2].

В окрестности разряда УПГ формируется импульсное электрическое поле высокой напряженности, что позволяет использовать этот разряд для воздействия поля на различные материалы и объекты. Кроме того, способ формирования разряда высоковольтными импульсами по-

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

С феноменологической точки зрения разряд УПГ представляет интерес в том смысле, что позволяет исследовать разрядные процессы при заданной положительной или отрицательной полярности прикладываемого импульсного потенциала, управлять амплитудой и частотой следования импульсов. В результате можно выделить эффекты воздействия этих параметров на разрядные процессы в различных условиях.

Исследования разряда УПГ проводились в воздухе при давлении р = 1 Тор. Разрядная трубка имела длину I = 2.8 м, внутренний радиус гт = 9.5 мм. На конце трубки был расположен покрытие-электрод (рис. 1) с длиной 1ПЭ = 2 см из мелкой металлической сетки с ячейкой 1х1 мм. На покрытие-электрод с генератора подавали высоковольтные импульсы потенциала отрицательной полярности амплитудой фпэ = 10 кВ.

В данной работе основное внимание уделялось поведению профиля свечения разряда УПГ по его длине. Как видим (рис. 1), разряд УПГ имеет в трубке коническую форму, сужающуюся к концу разряда. Причем часто длина разряда меньше длины трубки. Профили свечения разряда определялись по цифровым фотографиям пу-

Рис. 1. Разряд униполярного пробоя газа. Покрытие-электрод в виде сетки расположен слева.

Давление р = 1 Тор, воздух

Рис. 2. Профили свечения разряда УПГ на разных расстояниях от ПЭ:

1 - 10 см, 2 - 110 см, 3 - 210 см. На графике профили нормированы по амплитуде

тем усреднения по длине 2 см. Их вид менялся по длине разрядной трубки (рис. 2) - полуширина профиля уменьшалась по мере удаления от покрытие-электрода.

Описание модели

В качестве модели разряда УПГ рассматривается плазма с плотностью зарядов пе в длинной цилиндрической трубке, в которой заряженные частицы появляются с частотой ионизации V., распределяются по сечению трубки под действием амбиполярной диффузии с коэффициентом Da, и гибнут на стенках. Коэффициент Da предполагается постоянным, а частота ионизации V. = v(r, ^ может зависеть от радиальной координаты г по сечению трубки и времени t.

Уравнение баланса заряженных частиц в плазме разряда для некоторого сечения трубки запишем в виде ([3, с. 62]):

д пе _ д2 п д t

С г

В данном уравнении изменение концентрации электронов пе по времени t, заданное в левой части, определяется процессами диффузионного ухода с коэффициентом Da и их гибелью на стенках трубки (первое слагаемое справа) и рождением зарядов в объеме трубки с частотой ионизации V. (второе слагаемое справа).

Коэффициент амбиполярной диффузии Da зависит от подвижности ионов ¡и+ и температуры электронов Те и для неизотермической плазмы газового разряда равен ([3, с. 60]):

о.-=,,к-±,

е

где к = 1.38-10-23Дж/К - постоянная Больцмана, е = 1.6-10-19 Кл - заряд электрона. Для характерных условий разряда в воздухе (давление р = 1 Тор, Те = 5 эВ) коэффициент амбиполярной диффузии составляет Da« 2-103 см2/с.

Частота ионизации V. определяется энергетическим спектром электронов и зависит от потенциала ионизации I и плотности молекул или атомов газа п. Для максвелловского распределения электронов по энергиям частота ионизации V. равна ([3, стр. 69]):

^=п^СД + 2кТе )ехр(

кТе

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

V =

- тепловая скорость электронов.

Для основных составляющих воздуха - азота и кислорода в широком диапазоне температур электронов коэффициенты С., потенциалы ионизации I и вычисленные частоты ионизации V. для давления р = 1 Тор и температуры электронов Те = 5 эВ представлены в таблице 1.

Введем характерные параметры. В качестве характерной длины примем диаметр трубки d0 = 2 см. Тогда безразмерная радиальная координата будет равна р = г / d0, а концентрация элек-

Таблица1

С, 1017 см2/эВ I, эВ V, с'1

N2 0.85 15.6 4.7107

02 0.68 12.2 6.5107

тронов и = п-(dJ3. За характерное время возьмем длительность импульса тд, которая для наших условий составляет примерно тд = 4Ш6 с. Тогда безразмерное время будет равно т = t / тд, а безразмерная частота ионизации q = н.тд. Безразмерный коэффициент амбиполярной диффузии запишется в виде D = £атд / .

Заменим все переменные уравнения баланса заряженных частиц в плазме разряда на безразмерные величины:

д

d,

V“ о

D ■ d2

V™ 0

и q

2 + —3 — д(т ■ т0) т0 д(р ■ d0) d,

0

т

сократив затем каждое слагаемое на множитель

1

d3 ■ т

ио 10

, получим уравнение баланса заряженных ча-

стиц в плазме разряда для безразмерных величин:

д и д т

- = £-

д 2и ~др2

Искомая функция и = и(р,т) зависит от радиальной р и временной т координат. На расчетном отрезке [0, I] значения функции ионизации q = q(р,т) зависят от радиальной р и временной т координат. Задавая эту функцию и ее поведение во времени, будем добиваться распределения искомой концентрации заряженных частиц и = и(р,т) на отрезке с С [0, I], соответствующем наблюдаемому в эксперименте. При моделировании можно использовать аппроксимацию на сетке в виде явной или неявной схемы Эйлера.

Явная схема Эйлера

Безразмерное уравнение баланса заряженных частиц имеет вид уравнения теплопроводности. Его будем решать методом конечных разностей [4, с. 228]. Для этого разобьем отрезок [0, I]

Рис. 3. Конфигурация узлов сетки для явной схемы Эйлера

на N частей с координатами точек разбиения h¡ = і-Ш, где і = д..^ и шагом разбиения Лh = Ш. Введем шаг по времени Лт. Количество шагов по времени обозначим Nt, а номер текущего шага

k = д...т.

Непрерывные искомую функцию и(Ъ,т) и функцию ионизации q(h,т) заменим на дискретные функции и,к = и(Ъ,тк) и q¡к = q(h.,т|), заданные в точках (Ъ, тк). Затем заменим временную и пространственную производные, входящие в уравнение на соответствующие конечные разности (рис. 3):

д и Ли,

д т Лт

иі,к+і - иі, Л т

д2 и Л(Лиік)

иі+1,к 2 ■ иі,к + иі-1,к

(Л Ъ)2 (Л Ъ)2

,к+1 - искомая функция в следующий

д Ъ2 Здесь и

момент времени через малый интервал Лт, а и.к -в предыдущий. В результате уравнение баланса заряженных частиц в конечных разностях будет иметь следующий вид:

иі,к+і - иі, Лт

(ЛЪ)2

для всех і = 1..^ - 1.

Выразим в данном уравнении и.к:

—

иі,к+і =—■{иі-і,к+им,к)+ [і - с + qi,k ■ Лт\

где С

2 ■ £ ■ Лт (ЛЪ)2

- коэффициент Куранта, опре-

деляющий устойчивость разностной схемы. Начальные условия: и, = 0.001, для всех i = 0...Ж Граничные условия на концах отрезка при /' = 0 и /' = N зададим в предположении гибели зарядов у стенок разрядной трубки:

и, = 0 и и.„ = 0 для всех k = 0... Nt.

0,к ^к ^

Полученная схема расчета позволяет вычислять последовательно с заданным шагом по времени Ат распределение плотности электронов по сечению разрядной трубки и.к для всех /' = 0..М.

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

и

к

много меньше ее значения. При этом чем меньше выбирается шаг Аh разбиения отрезка [0, I], тем меньше должно быть значение шага по времени Ат. Принято, что вычислительная схема будет устойчива, если коэффициент Куранта меньше единицы С < 1. Шаг Аh входит в знаменатель коэффициента Куранта в квадрате, поэтому при уменьшении Аh в 2 раза Ат необходимо уменьшать как минимум в 4 раза для обеспечения устойчивости вычислительной схемы.

Неявная схема

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

Чтобы построить неявную разностную схему для уравнения диффузии, для дискретизации пространственной производной будем брать значения сеточной функции с верхнего (неизвестного) слоя по времени (рис. 4). Таким образом, разностное уравнение для (,, к)-го узла будет отличаться от уравнения для явной схемы только индексами по временной координате в правой части:

и1+1,к+1 2 ' и1к+1 + и1-1,к+1

+ q¡,k • и,к

Ат (А к)2

для всех /' = 1..М - 1.

Если привести подобные слагаемые, то получится система уравнений, связывающая для каждого 1-го узла три неизвестных значения сеточной функции (в самом этом узле и в соседних с ним слева и справа узлах). Множители при неизвестных значениях сеточной функции получим подобно тому, как это было сделано для явной схемы:

С (

- ■\и1-1,к+1

= [1 + Ч,к •Ат]•

1 )-(1+с )•

для всех /' = 1..М - 1.

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

Метод прогонки

Для решения систем линейных алгебраических уравнений с ленточной (в простейшем случае - трехдиагональной) матрицей разработан специальный алгоритм, называемый методом прогонки [4, с. 121], который является не чем иным, как методом Гаусса, учитывающим структуру матрицы. Суть метода прогонки состоит в том, что, как и в методе Гаусса, с помощью очередной строки производят исключение (обнуление) поддиагональных элементов столбца. Однако, поскольку известно, что матрица имеет ненулевые элементы, расположенные вдоль диагонали на ленте шириной 3, процесс Гаусса не доводят до конца, а ограничиваются лишь обнулением первых поддиагональных элементов.

Система уравнений для трехдиагональной матрицы имеет вид:

0

А

g1 g 2

г—1Л+1 /,¿+1 1+1,£+1

Лк

ЛЬ

1,к

Рис. 4. Конфигурация узлов сетки для неявной схемы

Прямым ходом прогонки мы обнуляем расположенные под главной диагональю элементы а1№ а21, аи1, ..., апп_1 и получаем двухдиагональную матрицу:

0

~2

и

1.к

ии+1 - ии

а

аю ап ап

и

а 21 а 22 а2з

и

2

0

а11 а12

и

0

а22 0^23

и

2

0

V... у

Элементы, помеченные тильдой, преобразуются в результате прямого хода. Формулы преобразования имеют вид:

аи = аи — ai-uвi, 8і = 8і — 8-і

и рассчитанное распределение плотности зарядов вдоль трубки и,к принимается за исходное

в.=-

І —І.І —I

І = 1,... П.

С помощью обратного хода по следующим формулам получаем решение системы линейных уравнений:

и = ~ и = ~ - аа+1и+1 ■ V 1 П

ип = , и-----~----, I = N - 1, ..., 0.

аа

Алгоритм вычислений

При организации вычислений необходимо задать начальное распределение зарядов по радиусу трубки. Мы предполагаем, что в начале расчета заряды равномерно распределены по всему диаметру трубки с некоторой начальной плотностью и,д, которая достаточно мала. Кроме того, требуется задать вид функции ионизации q¡k как в начальный момент времени, так и в течение всего времени численного эксперимента. Ее вид и поведение зависят от модели процесса ионизации. В случае конечных размеров области ионизации она может быть прямоугольного, треугольного вида или другой формы. Область ионизации может быть или постоянной или меняться в течение времени. Эти параметры выбираются так, чтобы результат численного моделирования соответствовал наблюдаемому в эксперименте. В результате можно сделать выводы о характере процессов ионизации в реальности.

Собственно, вычислительный процесс заключается в организации цикла, на каждом шаге которого время увеличивается на некоторую малую величину Ат и в следующий момент времени вычисляется новое распределение плотности зарядов по сечению трубки и,к+1, получающееся в результате совместного действия процессов рождения, гибели и диффузии зарядов в плазме разряда униполярного пробоя газа в длинной цилиндрической трубке исходя из предыдущего распределения и,к. Вклад этих процессов учитывается с помощью задания соответствующих коэффициентов: амбиполярной диффузии Б, а также формы и поведения функции ионизации q. После прохождения цикла вычислений учитыва-

ются граничные условия: и

0 и и.

0

и., = и.,^, для всех і = 0...Ы.

І^ І^+1

Обсуждение результатов

Известно [5], что поведение свечения разряда отражает поведение его заряженной компоненты. Поэтому расчет плотности зарядов в разряде можно при моделировании сравнивать с распределением интегрального свечения. Задачей данного исследования было выяснить, какие процессы в первую очередь определяют наблюдаемое на эксперименте поведение интегрального свечения разряда как по сечению разрядной трубки, так и в течение импульса. Фотографирование разряда цифровой камерой позволило получить радиальное распределение интегрального свечения, усредненное по времени импульса (рис. 2).

Выбор поведения во времени сеточной функции ионизации q¡k определяется экспериментом. В нашем случае разряд возбуждался импульсами потенциала отрицательной полярности. При этом даже при достаточно большой длительности импульсов (до 20 мкс) длительность ввода энергии в разряд, которую можно оценить по длительности фронта импульса тока, составляла порядка

0.5 мкс. В течение этого времени энергия, вводимая в разряд, сначала нарастает, затем уменьшается. Такое поведение можно смоделировать функцией Гаусса. При моделировании используется последовательный набор профилей функции ионизации q¡k через временной интервал Ат. Для сравнения с экспериментом рассчитанные для разных моментов времени профили были просуммированы и нормированы на 1.

В качестве радиального распределения сеточной функции ионизации q¡k можно было использовать прямоугольное, треугольное и обратное. Прямоугольный профиль соответствует одинаковой скорости ионизации по сечению трубки, треугольный - преимущественной ионизации на оси трубки, а обратный - у ее стенок. Сравнение результатов моделирования и измеренных экспериментально профилей свечения разряда УПГ показало, что рассчитанный радиальный профиль свечения соответствует наблюдаемому в эксперименте. Причем вблизи покрытие-электрода ионизация происходит в основном у стенок трубки (рис. 5а), в середине разряда - наблюдается согласие при условии равномерной ионизации по всему сечению трубки (рис. 5б), а в конце раз-

І.І—1

а) б) в)

Рис. 5. Сравнение рассчитанных и измеренных (+) радиальных профилей разряда УПГ на разных расстояниях от ПЭ: а) 10 см (пристеночная ионизация), б) 110 см (равномерная ионизация), в) 210 см (осевая ионизация). Жирная линия - распределение функции ионизации. Тонкими линиями показан временной ход профиля свечения в течение импульса

ряда - только при условии преимущественной ионизации на оси трубки (рис. 5в).

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

Библиографический список

1. Герасимов И.В. Излучательные свойства разряда униполярного пробоя газа // Журнал технической физики. - 1994. - Т. 65. - С. 30-35.

2. Савинов В.П., Сухов А.К., Копейкина Т.П. Изменение функции распределения электронов по энергиям в плазме разряда униполярного пробоя газа // Вестник Московского университета. Сер. 3. Физика. Астрономия. - 2009. - Т. 50. - №9 3. -С. 56-59.

3. Райзер Ю.П. Физика газового разряда: Учеб. руководство: Для вузов. - 2-е изд. перераб. и доп. - М., 1992. - 536 с.

4. Турчак Л.И., Плотников П.В. Основы численных методов: Учеб. пособие. - 2-е изд. перер. и доп.- М., 2003. - 304 с.

5. Королев Ю.Д., Месяц Г.А. Физика импульсного пробоя газов. - М., 1991. - 224 с.