Двумерное моделирование формирования тлеющего разряда в аргоне. Часть 1. Модель и расчетная схема Текст научной статьи по специальности «Физика»

УДК 533.9

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

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

[email protected]

ДВУМЕРНОЕ МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ ТЛЕЮЩЕГО РАЗРЯДА В АРГОНЕ. Часть 1. Модель и расчетная схема

Для компьютерного моделирования двумерного распределения параметров плазмы в тлеющем разряде в аргоне использовалась модель плазмы газового разряда в двумерных декартовых координатах. Система уравнений двух-жидкостной гидродинамики содержала два уравнения непрерывности для электронной и ионной компонент плазмы и уравнение электростатики в виде уравнения Пуассона для вектора напряженности электрического поля E.

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

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

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

Газовые разряды широко применяются на практике в плазмохимических реакторах, резке металлов, мощных источниках света, для накачки лазеров, в ионных двигателях космических аппаратов и т.п. Однако степень изученности всего многообразия газовых разрядов не является достаточной. Например, слабо изучены униполярные формы разряда при низком давлении с асимметричной конфигурацией электродов [1-5]. Важным инструментом исследования процессов в газовых разрядах является их моделирование. Хотя модель использует некоторые упрощения реальных процессов, моделирование позволяет выявить определенные особенности и эффекты, проявляющиеся в реальных разрядах [6-7].

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

Так, в работе [10] двумерное моделирование использовали для расчета тлеющего разряда постоянного тока в условиях существования нескольких токопроводящих столбов. Посредством численных экспериментов исследован процесс объединения двух столбов непрерывного тлеющего разряда. В работе [11] трехмерная модель применялась для расчета пространственной конфигурации химически активного тлеющего разряда в воздухе с бес-

конечными плоскопараллельными электродами. Получены пространственные поля распределения электронной температуры и заселённостей колебательных уровней азота. Двумерное моделирование тлеющего разряда в аргоне сложной формы проведено в работе [12]. Расчеты показали, что схождение решения системы исходных уравнений к стационарному состоянию очень чувствительно к первоначальной оценке для значений неизвестных параметров. Поэтому для обеспечения сходимости подбирали начальные условия. Показано хорошее согласие расчетов с экспериментальными результатами.

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

Задачей данной работы, которая является логическим продолжением исследований [6-7], ставится двумерное моделирование поведения пространственного распределения параметров тлеющего разряда в аргоне в ходе его формирования.

© Сухов А.К., 2014

Вестник КГУ им. H.A. Некрасова «й- № 6, 2014

19

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

Модель разряда задана системой дифференциальных уравнений двухжидкостной гидродинамики в двумерной декартовой геометрии, описывающих поведение параметров плазмы с учётом граничных условий в приэлектродных слоях и на стенках [11-13].

Разряд описывали системой векторных уравнений непрерывности для электронной и ионной компонент и уравнением Пуассона для вектора напряженности электрического поля Е:

+ Шу Г = а\Т \-Bnn ,

л, е е / е р >

Т е =--е^еЕ - Ъягайпе, д- - 1-1

— + Тр =а\Те\ в—е— р , (1),

Тр = -рМрЕ - -р ,

Шу Е = 4ле(-р - -е).

Е = - grad ф

где -е и -р - концентрации, Т е и Тр - потоки, це и цр -подвижности, Б е и Бр - коэффициенты диффузии, соответственно, электронов и ионов, Е - напряженность электрического поля, а - первый ионизационный коэффициент Таунсенда, в - коэффициент объемной рекомбинации, ф - потенциал, е - заряд электрона. Б е = це-Те и Бр = ЦрТр - коэффициенты диффузии, соответственно, электронов и ионов.

В двумерном случае декартовых координат система (1) принимала следующий вид:

^ + ^ + ^ = а\Т I-Р-п , дt дх ду

Тех -еЕх + Те

(

Гу =

J7 , Т дПе

n E + T —1

1 у 1 дУ

дп дГ дГ

—— + ^ + ^ = а\Г I -впепр д дх ду

(2).

(

Г Рх = V—

Г—У = Vp

дп р ^

ПрЕх - Tp^~

v дх /

( дпр}

п—Еу- T— -у

дпе

3t

+ Ve

д ф дф дп д пе

п —— +—---1— T --

1 дх2 дх дх e дх2

д2ф дф дпе д пе + п —— +—---1— 1 --

1 ду2 ду ду 1 ду2

= а\Ге\ -Рп1пр>

-V—

дф+дфдп^ + т р дх2 дх дх р

д\ дх2

(3);

д2ф дф дп— д2п— + пр—т +——- + T —f-ду ду ду ду

= а\Г J -вп1п—>

д 2ф д 2ф . ,

здесь модуль потока электронов равен:

Г =Лп ф-T

х

V дх дх

дф

дп„

Г1у=vV - T-t

Г 1 =V Г 2х + Г 2у.

Начальные условия

Концентрации электронов и ионов в начальный момент времени во всей расчетной области считали неизменными и равными фоновой концентрации п0:

пе(x, y, 0) = np(x, y, 0) = n0.

Градиент потенциала в направлении оси 0х в начальный момент времени считали постоянным, в направлении оси 0у - равным нулю:

дф(х, у, 0) = const; — (х, у, 0) = 0.

дх ду

Граничные условия

Потенциал на катоде 1 считали равным нулю, а на аноде 2 - приложенному напряжению V, здесь d - длина разрядного промежутка:

Ф(0, y, t) = 0, q>(d, y, t)= V.

Градиент потенциала в центре и у стенок считали нулевым, а у катода 1 и анода 2 - постоянным, здесь r0 - размер разрядного промежутка по оси 0у:

дЕ дЕу

—Г + = 4Ж (-р - -е ^ х у

Ех =-Ф, Е = -ф. дх у ду

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

Рис. 1. Расчетная область (К - катод, А - анод, 0х - ось разряда)

+

+

0 = о, p, г0, t) = 0,

ду ду

дР (x, 0, t) = const; dpp (x, d, t) = const.

dx dx

Граничные условия на электродах задавали исходя из того, что электроны испускаются катодом под действием падающего потока ионов Г (0, y, t) с коэффициентом ион-электронной эмиссии у, а анод ионы не испускает:

Гех(0, y, t) = -уГ/0, y, t);

rjd, y, t) = 0.

В этом случае краевые условия для плотностей зарядов на катоде 1 и аноде 2 записывались в виде:

n(0, y, t) = у np(0, y, t)-^p/pe ;

n(d, y, t) = n0 .

На оси трубки 3 плотность зарядов считалась непрерывной, а на стенке 4 - фоновой:

dn dn

(x,0, t) = 0;-^- (x,0, t) = 0;

ду ду

ne(x, rff t) = no; np(x, r„ t) = no .

Для решения заданной системы уравнений использовали численные методы. В расчетной пространственно-временной области 0 < x < d, 0 < y < r0, 0 < t < z nk вводили конечно-разностную сетку: mht = {x. = ihi, i = 0 ... ni; y. = j hj, j = 0 ... nj; t* = кт, k = 0 ... nk} с пространственными шагами по координатам hi = d/ni и hj = r</nj и шагом по времени т. Здесь , i = 0 ... ni; j = 0 ... nj - индексы по координатам, а k = 0 . nk - индекс по времени.

Рассматривали два временных слоя: нижний t* = кт, на котором распределение искомых функций известно и верхний временной слой t*+1 = (к+1)т, на котором распределение искомых функций подлежит определению. Непрерывные функции потенциала, потока, концентраций и поля заменяли сеточными: ф(х, y, t*) = m. *, Г (x, y, t*) = Г .*, Г (x, y, t*) = Г *

' j i,j ex i^j exi,j 7 eyx r ^ j / ey

i,j 7 exv r ^ j

Г (х,у, Л = Г .к, Г (х , у, Л = Г .к, п (х ,у, Л = п ..к,

рх Г у у у рх1,] РУ Г I РУ1'У е г ^ у / е1,]

п (х, у, 1к) = п .к, Е (х, у, ?) = Е.к,

р Г I Р1,] х I у у / XI,] 7

Е (х , у, ¡к) = Еуук. Частные производные, входящие в уравнения, заменяли конечно-разностными соотношениями. Для временных производных использовали правые, а для пространственных - центральные разности:

Полученная система уравнений в конечных разностях записывалась в виде:

*+1 * n . . - n . .

ei,j ei,j

+ Ue

Р+ц - 2<p*j + pj-1J hi 2

* PP.+i - 2pp. +Ppj-i + nij -Z/T-— +

hj2

+ (P*i,j -PP-l.j )-(ne*+1,j - ne*-1,j ) 4hi2

(ppj+1 -Ppj-1 )'(ne*j+1 - ne{j-1 )

4hj2

* ~ * *

^ , - 2^ . + n . , .

T ei+1, j ei, j ei-1, j

- T

hi2

* - * * ne - 2ne + ne

hj2

* t-I * о * *

= a. .U . - en . n

l, j| ei,j\ ^ ei, j pi

*+i * n - n

Pi,j Pi,j

* Pp+i,j - 2ppj + P*-i,j n ---i1-— +

pi,j hi2

* Ppj+i - 2ppj +ppj-i + nP*j-~hf2-+

(Pi+i,j -Pp-ij)-(nph,j -np*-hj)

+----J— _

4hi2

-i HnpLi- nP*,j-i)

(ppj+i-p, j-i )-(nPi,j+i- npiJ-l,

+-г-; ^ -— +

4hj2

* ~ * *

n - 2n + n

+ T pi+1,j pi,j pi-1,j

+ p hi2

* n, * It

n - 2n + n

Pi.j+1 Pi,j Pi.j-1

+T

hj2

* Г"» * n * *

= «i,j|Гei,j\ -PneUinpUj

pp+hj - 2ppj +pP-hj

hi 2

ppj+i - 2p*j +ppj-i +----jr-— = 4-jieu. .. .. .

hj2 ei,j Pij

i = 1...ni -1, j = 1...nj -1, * = 0...n* -1.

= 4ne(n . . - n ),

v ei, j Pi, l''

(4).

После ее упрощения получили систему следующего вида:

= ne*j +T[dnaaj - faph -

ei, j ei,j

- и (dn* . .+ dn* )l

n„*+1 = n* + r\dna*j - dnB*i +

i, j i, j

Up (dn*px.j + dnPyiJ )j

Pp+1,j - 2ppj + Pp-1,j + Ppj+1 - 2ppj + ppj-1

(5),

hj2

hi2

= 4ne(ne* , j - npI j ),

i = 1...ni -1, j = 1...nj -1, * = 0...n* -1.

к

+

n

i, j

X

z

+

где:

dnak, = ak, Г к

'j Я

к

^Pij =Рп„,1при, ¿п*лj =[пе*j(ф+1,j -2pPj j)+

+ 0,25фj -фф-1,j)-(пек+1,j - п.;-,,j)-

Vi.i+1,j - пи- Te (пк,, , - 2пк, + пк, ,.)]/hi2

1 i -1, j

j=k jфj+i - ф+ppj-,)+

+ 0,25(ppj+i -i- nj)-- T (п* .+, - 2п1к. + п1к . , )|/hj2,

j+1 ei,j ег,j-1/J J 5

j =[п—кj Фj - Pj + Pi,j )+ +0,25(Ф+1,^ ф^)-(п—к+1,j - п—к-i,^)

(6).

+ Tp(ппк, - 2п„к + пк, )|/hii

p\ Pi+i,j Pi,j —i-i,j

Kj =k", jфj+1 - 2P + pPj-1)+

+0,25(ФФ,j+i Ф у(п—к м - п;1м)

,(п„к - 2прк + пк )/ hj2,

Л pi,j+1 pi,j pi,j-l). J '

+ T

-Г к Ve \ к ( к к \

j = К j (Ф<+1, j -Ф-1, j )-

кк -T1 K+i,j -п

к

I,i )1

к V1 к к к Feyj = 2hj ^j((P',j+i)-

- T (fC,j+i - п1к,j-iJ

(7).

(

ai , = 28,5 • p • exp

- 25,9

Ек

здесь модуль электрического поля:

E =-*i, j

к

Pi+i,j ■

к

Pi-i, j

2hi

E к =

Еуи j

к

Pi,j+i ■

к

■Pu j-i

2hj

(8)

(9).

\Ек I = Ек2 + Ек2

ГМ \n'xi,Jn'Уi, j

Для вычисления модуля потока электронов использовали соотношения:

\гк 1= Гк 2 + Гк 2

| 1i, j\ \ ^i, j 1уi, j'

Первый ионизационный коэффициент Таун-сенда aj, зависящий от модуля напряженности электрического поля |Еук| и давления аргона р (при T = const), аппроксимировали следующим выражением [i3], в котором коэффициенты подобраны по данным работы [i4]:

Рис. 2. Вычислительная блок-схема расчетов

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

Метод верхней релаксации

Для расчета распределения потенциала в разрядном промежутке использовали метод верхней релаксации [15], который в данных условиях давал наибольшее быстродействие. Двумерное уравнение Пуассона в декартовых координатах для конечных разностей имело вид: Ф* - 2ф*. 1

hi2

Ppj+i - 2ррj +ppj-i

hj2

= 4яв(п к. - п к ),

v j pi, j''

(i0).

i = i...m - i, j = i...Hj - i,

Коэффициент объемной рекомбинации в = 6,510-7 см3/с, коэффициент ион-электронной эмиссии у = 0,07. Подвижности электронов

Ле = 3,3105/р (см2Тор)/(Вс) и ионов ¡л = 2,0103/р (см2Тор)/(Вс).

Блок-схема программы

Алгоритм проведения вычислений представлен на блок-схеме (см.: рис. 2). В программе предусмотрено сохранение промежуточных результатов

На границах расчетной области задавали следующие краевые условия. Потенциалы на электродах полагали равными 0 и V:

фк0,1 = 0 фкп,,1 = ^ } = 0~ГУ.

Градиент потенциала у стенок считали нулевым:

Фки,

фк ■,,

' i,nj-P

i = 0...ni.

Перепишем уравнение Пуассона в виде:

Рис. 3. График зависимости времени расчетов от параметра со. Параметры расчета: т = 600, щ = 100

1

1

1

a V+1 =-Lr^,1 ■ +-W+1, + ■ + ■ , -

hi hj 2^',j-1 hi 2^,+hJ h2%-j+1 к к , 1 1

1

hj2

- 4n(nei . - npij), ^ a = + j

здесь значения ф и ф встречаются в расчете раньше, поэтому берутся из следующего (к + 1) -го слоя, а ф и ф - позже и берутся из текущего к-го слоя.

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

к+1 ® ,„к+1 . ® ,„к+1 . ^ ,„к , ^ ,„к = Щ2+ ] ^ + + ] ^ +

+ (1- о) ■a -vlj -ю 4же(пек. - npk,),

(11).

При программной реализации данной формулы применяли внешний цикл по j от 1 до nj - 1 и внутренний по i от 1 до ni - 1. Параметр те[1, 2) выбирался так, чтобы время расчетов было минимальным. Чем больше было количество разбиений по координатам, тем сильнее оптимальное значение т сдвигалось к его верхней границе.

При соответствующем подборе ускоряющего параметра т можно было сократить время расчета в сотни раз. Для заданных параметров модели и разбиении сетки оптимальное значение было близко к т = 1,99.

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

1. Герасимов И.В., Сухов А.К. Исследование распространения разряда униполярного пробоя газа // Вестник Костромского государственного университета им. Н.А. Некрасова. Серия: Педагогика. Психология. Социальная работа. Ювеноло-гия. Социокинетика. - 2005. - № 11. - C. 5-10.

2. Герасимов А.И., Герасимов И.В., Сухов А.К., Якунина Л.В. Исследование спектра излучения разряда униполярного пробоя газа в воздухе при изменении давления // Вестник Костромского государственного университета им. Н.А. Некрасова. Серия: Педагогика. Психология. Социальная работа. Юве-нология. Социокинетика. - 2005. - № 3. - С. 10-14.

3. Герасимов И.В., Сухов А.К., Копейкина Т.П. Влияние частоты следования импульсов высоко-

вольтного потенциала на параметры разряда униполярного пробоя газа // Вестник Костромского государственного университета им. Н.А. Некрасова. - 2009. - № 3. - C. 30-34.

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

5. Сухов А.К. Две формы существования разряда униполярного пробоя газа // Вестник Московского государственного университета М.В. Ломоносова. Серия 3: Физика. Астрономия. - 2013. - № 2. - С. 56-60.

6. Сухов А.К. Моделирование радиального профиля свечения по длине разряда униполярного пробоя газа // Вестник Костромского государственного университета им. Н.А. Некрасова. - 2010. -№ 2. - С. 12-17.

7. Сухов А.К. Моделирование импульсного разряда в азоте // Вестник Костромского государственного университета им. Н.А. Некрасова. - 2013. -Т. 19. - № 2. - С. 15-24.

8. Райзер Ю. П. Физика газового разряда: науч. изд. - 3е изд., перераб. и доп. - Долгопрудный: Изд. дом «Интеллект», 2009. - 736 с.

9. Суржиков С.Т. Физическая механика газовых разрядов. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2006. - 639 с.

10. Суржиков С.Т. Численный анализ структуры двух типов тлеющих разрядов. // Физико-химическая кинетика в газовой динамике. - 2008. -Т. 7. - [Электронный ресурс]. - Режим доступа: http://chemphys.edu.ru/article/112/.

11. Петрусев А.С. Трехмерная численная модель для химически активного тлеющего разряда в воздухе // Труды МФТИ. - 2009. - Т. 1. - № 1. - С. 90-93.

12. Rafatov I.R., Akbar D., Bilikmen S. Modelling of non-uniform DC driven glow discharge in argon gas // Physics Letters A. - 2007. - Vol. 367. - P. 114-119.

13. Ward A.L. Calculations of cathode-fall characteristics // J. Appl. Phys. - 1962. - Vol. 33. - P.2789-2794.

14. Smirnov B.M. Reference Data on Atomic Physics and Atomic Processes. - NY.: Springer, 2008. - 173 p.

15. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978. -592 с.