CC BY
27
2017, Т. 159, кн. 4 С. 444-457
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
УДК 533.9.01
МОДИФИКАЦИЯ МЕТОДА ШАРФЕТТЕРА-ГУММЕЛЯ ДЛЯ НАХОЖДЕНИЯ ПОТОКА ЗАРЯЖЕННЫХ ЧАСТИЦ ПРИ МОДЕЛИРОВАНИИ ВЫСОКОЧАСТОТНОГО ЕМКОСТНОГО РАЗРЯДА
В.С. Желтухин1, М.С. Фадеева2, В.Ю. Чебакова2
1 Казанский национальный исследовательским технологический университет,
г. Казань, 420015, Россия 2Казанский (Приволжским) федеральный университет, г. Казань, 420008, Россия
Аннотация
В работе представлен численный метод решения нестационарного уравнения конвекции-диффузии для нахождения концентрации заряженных частиц, одновременно позволяющий рассчитывать значения плотности потока заряженных частиц при наличии больших градиентов как плотности заряженных частиц, так и коэффициентов уравнения. Метод получен путем модификации алгоритма Шарфеттера-Гуммеля. Сначала вычисляется концентрация заряженных частиц с помощью неявной разностной аппроксимации, построенной интегро-интерполяционным методом. Затем вычисляется плотность потока заряженных частиц. Предлагаемая модификация позволяет проводить расчеты при более слабых, чем условие Куранта, ограничениях на выбор шага временной сетки. Кроме того, допускается непостоянство коэффициента перед градиентом в выражении для плотности потока, а также наличие источникового члена.
Ключевые слова: математическое моделирование, высокочастотный емкостной разряд, модифицированный метод Шарфеттера-Гуммеля, уравнение конвекции-диффузии, плотность потока заряженных частиц
Введение
Низкотемпературная плазма широко используется при получении материалов с улучшенными свойствами [1—8]. Важную роль при этом играет применение высокочастотных емкостных (ВЧЕ) разрядов [9-13]. При моделировании ВЧЕ-разрядов для нахождения энергии, передаваемой электрическим полем электронам, возникает необходимость расчета плотности электронного потока при больших градиентах концентрации электронов, а также коэффициентов уравнения, описывающего баланс заряженных частиц. Кроме того, одной из основных характеристик каждого конкретного плазмотрона является вольтамперная характеристика (ВАХ) -зависимость между полным током через газоразрядный промежуток и напряжением на нем, построенная с помощью либо натурных экспериментов, либо численных расчетов. Полный ток = де(0+ — 0е + О2+) + еодЕ/дЬ, проходящий через газоразрядный промежуток, представляет собой сумму тока проводимости де(0+ — 0е + О2+) и тока смещения е^дЕ/дЬ. Здесь де - заряд электрона, 0е -плотность электронного потока, 0+ - плотность потока атомарных ионов, О2+ -плотность потока молекулярных ионов, ео - электрическая постоянная, Е - напряженность электрического поля [14-19]. Следовательно, для расчета полного тока необходимо также вычисление плотностей заряженных частиц, как электронов, так и ионов.
Вычисление плотности потока электронов Ое = —Педпе/дх — ¡лепеЕ требует применения численного дифференцирования. В случае сильных градиентов пе, Ве, Е, имеющих место в приэлектродных областях, возникает существенная погрешность в вычислении Ое. При решении уравнения баланса электронов эти трудности можно обойти с помощью метода потоковой прогонки [22], однако он неприменим при использовании направленных разностей, необходимых для вычисления Ое. В настоящей работе предложен численный метод решения, позволяющий вычислять плотность потока заряженных частиц при наличии сильных градиентов пе, Пе, Е. Метод основан на модификации явного двухшагового метода Шарфеттера-Гуммеля, изложенного в [20, 21]. Этот метод имеет весьма обременительные ограничения на шаг по времени. В предложенной нами модификации изменен порядок действий. На первом шаге вычисляется концентрация заряженных частиц с помощью неявной разностной аппроксимации, построенной интегро-интерполяционным методом с применением направленных разностей, что позволяет учитывать только физические соображений при выборе временного шага (например, время жизни метастабильных атомов). На втором шаге непосредственно вычисляется плотность потока заряженных частиц с учетом переменной подвижности заряженных частиц.
Отметим, что существует ряд методов, основанных на построении модифицированного функционала Лагранжа и позволяющих производить одновременное вычисление как самого решения, так и его градиента (см., например, [23-26]). Но данные методы неприменимы в рассматриваемом нами случае, поскольку они предназначены для решения в основном стационарных задач.
1. Постановка задачи
Рассмотрим ВЧЕ-разряд между двумя плоско-параллельными симметричными электродами, один из которых заземлен, а другой соединен с ВЧЕ-генератором, расстояние между электродами значительно меньше размеров самих электродов. Обозначим через Ь расстояние между электродами, считаем, что заземленный электрод находится в точке х = 0, нагруженный - в точке х = Ь (ось Ох направлена перпендикулярно поверхности электродов). Математическая модель ВЧЕ-разряда в нелокальном приближении при пониженных давлениях в аргоне [27, 28] содержит в том числе уравнение баланса электронов
дп д ( дп \
+ дх \ дх — п^еЕ ) = ЕпК + Е2пт + Язптпе — Я4пеп+ —
— Еъп+ п2е + Яюп2т — Нт+пе + + Нцптп2*, 0 < х < Ь, г > 0, (1)
где пе, п+, пт, Н2*, п2+ - концентрации электронов, атомарных положительно заряженных ионов, метастабильных атомов аргона, димеров и молекулярных положительно заряженных ионов аргона, соответственно, Ое = —Ве дпе/дх — — пе^еЕ - плотность потока электронов, ¡ле - подвижность электронов, Ве - коэффициент диффузии электронов. Здесь и далее Н, г = !,..., 5,10,13,16,17, -коэффициенты скоростей плазмохимических реакций Аг * +Аг * ^ Аг + +е (хемо-ионизация); 2Аг + Аг + ^Аг + +Аг (конверсия ионов); Аг + +е ^ Аг + Аг + +е; Аг + + е ^ Аг + Аг; Аг * +Аг * ^ е + 2Аг + Аг + ; Аг 2 +Аг * ^ е + Аг + Аг + ; е + Аг * ^ 2Аг* +е. Значения коэффициентов этих плазмохимических реакций приведены в [27, 28].
Граничные условия при х = 0 выбираются в виде
(Се(0,г) = (—пеУе — ЦепеЕ)\(0+^) , Е > 0,
1 (2)
[Се(0,£)= (—пРУе — 7(^+п+Е + 1Л2+п2+Е))\(0+^), Е< 0;
а при х = Ь - в виде
(0е(Ь,г) = (ие,уе — /лепеЕ)\(Ь_ ^ , Е < 0, (3)
[0е(Ь,г)= (ПеУе — 7(^+П+Е + Ц2 + П2+Е))\(Ь-,4) , Е > 0.
Здесь уе = у^8кТе/(пте)/4 - средняя тепловая скорость электронов, 7 - коэффициент вторичной электронной эмиссии с электродов, Те - температура электронов, те - масса электрона, Ц2+ - подвижность молекулярных ионов.
Во многих случаях величина пе изменяется очень сильно в расчетной области вблизи электродов. Поэтому вычисление плотности потоков заряженных частиц по приближенно найденным значениям пе , как правило, приводит к существенным погрешностям. В настоящей работе предлагается модификация метода Шар-феттера - Гуммеля [21], позволяющая проводить вычисление плотности потоков заряженных частиц с приемлемой точностью.
Отметим также, что если вычисление концентрации пе проводится с помощью явной разностной схемы, как в методе Шарфеттера-Гуммеля, то для обеспечения устойчивости схемы необходимо выполнение условия Куранта. При наличии больших полей Е данное условие требует использования очень маленького временного шага, что существенно увеличивает время расчетов. Кроме того, наличие источ-никового члена, осциллирующего с частотой генератора поля и отсутствующего в [21], приводит к дополнительной неустойчивости вычислений. Использование же неявной разностной схемы позволяет избежать указанных трудностей.
2. Вычисление концентрации заряженных частиц
Разностную аппроксимацию для уравнения баланса (конвекции-диффузии) электронов (1) будем строить интегро-интерполяционным методом [29] с применением направленных разностей [30]. Введем на отрезке [0, Ь] пространственные равномерные сетки йн = {х1 = 1Н, I = 0,1,... ,М }, йн = {х1 = 1Н, I = 1, 2,... ,М — — 1}, Н = Ь/М - пространственный шаг, а также временную сетку йТ = = = вт, в = 0,1,...}, т = Т/Ьт - шаг по времени, Ьт - количество временных шагов на периоде Т осцилляции электрического поля. На сетках йн х йТ, йн х йТ определим сеточные функции (для которых сохраним те же обозначения, что и для дифференциальных функций). Значения сеточных функций р в точке х® = (х1,Ье) € йн х йТ будем обозначать через .
При построении разностных схем значения коэффициентов дифференциальных уравнений опускаются на нижний временной слой \. Проинтегрировав уравнение (1) по отрезку [х1 — Н/2, х\ + Н/2] и полагая Ь = , получим
Х1+Н/2 хг+Н/2
/дп Р д / д п \
+ дь\ — ^е + 1^еПе(—Е)) (1х = =
х\-Н/2 х\-Н/2
Х1 + Н/2
= ! КгПеМ + К2П2т + Язптпе — Я^ПеП+ — Я5п+п2е+
XI -Н/2
+ Я\оп2т — Я\зП2+Пе + Ягеп^ + Я17птп2* ¿х. (4)
Интеграл . заменим квадратурной формулой центральных прямоугольников, при этом производную по времени будем аппроксимировать разностью назад.
Далее,
Х1 + Н/2
■72 = ! — °е + Мепе(—Е)) ¿х
Х1-Н/2
= (°е д~х) — В д~х)
дх ) (х1-к/2,ьв) \ дх ) (Х1+к/2,ьв)
) — В дп) + (5)
/ (х1-к/2,гв) \ дх ) (х1+к/2,ьв)
Х1 + Н/2
_ ___ дпе
V е дх ) (х1-н/2,г3) Vе да
д^РпРХ—Е)
где J21 = J -д-^х. Первые два слагаемых в правой части последнего
Х1-Н/2
выражения аппроксимируем разностными отношениями / дп \
[Ве-^] , , « Ве,1±1/2пе,х(х1±1/2,Ьз).
\ дх ) (Х1±К/2,ЬВ) ' '
Для аппроксимации 321 предварительно преобразуем подынтегральное выражение, выделив при этом у функции Е положительную и отрицательную составляющие: Е = Е + + Е- , Е± = (Е ± \Е\)/2. Получим
+
д^пе(—Е) д(—Е) д^е'п-е д(—Е)
-д- = Меп^ —--+ ( — Е )—- = Мепе -д--+
х х х х
. ( 1?)+ дМепе , д(—Е) ( д Мепе , рЛ + ( — Е Г^- + Мепе-д- + (—Е) —К-. (6)
х х х
Тогда
Х1 + Н/2 Х1 + Н/2
Т± [ д(—Е)±1 , [ ( Т?)± дМепе , т± . т± ^21 = Ме пе -дх- + (—Е) дх = 3210 + 3211.
Х.1-К/2 Х.1-К/2
Интегралы приближаем следующим образом:
3±10 х (Репе)
Х 1+Н/2
I х
ХI ] дх
Х1-К/2
= (Ме пе) [(—Е)±(х1 + Н/2,13-1) — (—Е)±(х1 — Н/2,13-1)]. (7)
ХI
Х ь + Н/2
/д м п
(—Е)± —^ е ¿х приближаем с использованием разно-
Х1-Н/2
стей против потока, при этом считаем, что направление движения потока электронов обратно направлению поля Е , причем за положительное направление поля Е выбираем положительное направление оси абсцисс. Следовательно, в случае отрицательного направления движения потока электронов производную д(ре пе)/дх заменяем разностным отношением
(Ме пе) х(хь ,Ьв) = [МеЫ^в-^пе^Ьэ) — Ме(х1-1,Ьз-1)пе(х1-1,Ь3)] /Н. (8) В случае положительного направления движения потока электронов полагаем
д(Мепе)/дх X (Мепе) х(х1 + 1,Ьв) = [(М еп е(х 1+1, Ь — (Мепе)(хь,Ьв)] /Н.
Таким образом, интеграл ^ц заменяется квадратурной формулой левых прямоугольников, а интеграл ..—ц - правых прямоугольников, то есть
■ л ~ ( — Е) + (х1 — Н/2) [pe(xl,ts-l)ne(xl,ts) — ^le(xl-1,ts-l)ne(xl-1,ts)] , (9)
11 ~ (—Е)-(Х1 + Н/2) [^e(xl + l,ts-l )ne(xl + l,ts) — ^еЫ ,ts-l)ne(xl . (10)
Наконец, интеграл ■, как и интеграл ■1, заменим квадратурной формулой центральных прямоугольников
■ « Н
Яl(xl,ts-l)ne(xl ,ts)N (xl,ts-l) + ЩХ1 ,ts-l)nm(xl ,ts-l) + + ЯЗ(Х1 ,ts-l)nm(xl,ts-l)ne(xl ,ts) — ЯA(xhts-l)ne(Xl ,ts)n+(xl,ts-l) — — Я5(xl,ts-l)n+ (xl,ts-l)n2e(xl,ts) + Яlо(xl,ts-l)n2m(xl,ts-l) — — Я1З(Х1 ,ts-l)n2+(xl,ts-l)ne (xl,ts) + Я1ех ^-1)П2*(Х1 ^-1)+
+ Яlr(Xl,ts-l)nm(Xl,ts-l)n2*(Xl,ts-l) . (11)
При аппроксимации квадратичной нелинейности в правой части уравнения используется линеаризация по Ньютону вида
п2е (Xl,ts) « n2e(Xl,ts-l) +2пех ,ts-l) [ne(Xl,ts) — ^^ ^-1)] .
В результате разностная схема для уравнения конвекции-диффузии электронного газа имеет вид
п'е, l пе- , /п s — 1 пе, ^1 пе, l г^-1 , l пе, — 1
_Н- Ds-1 е'1,+1_еА _ ns-
т Н е^+1/2 Н Бе,— 1/2
( I7^^ + ?S — 1 S —1 | / ТГ^Л — ?S — 1 ^ — 1 | / T7^^ + JS—1 ^ — 1 s
— ( — Е) — 1/2 ^e,l-1ne,l-1 + (—Е)1 + 1/2 ^е^ + ^е^ + (—Е)1+1/2 ^ пе^ —
— (—Е)—-—1п^ = н[Я1—1п^ N-1 + Я2—'(п^Г1 )2+ + Я3— 'К—К! — Я- 4,1 п+— — Щ—1п+—11 [«— ')2 + 2п8е— 1(пе^ — пе—+
+ Я{—(п—)2 — Я—пЪ+пе^ + Щ-т-!)2 + . (12)
Для аппроксимации граничных условий (2) при х = 0 в случае симметричных электродов проинтегрируем уравнение (1) по отрезку [Н/2, х1 + Н/2], тогда получим
_п^-1 / Л г
Н + О-/ — 0—/2) = Н[я1-1пег 1Щ-1 + Я2-1(п—1)2 + + Я^-^пГ—^^ < 11 — ЯА—11п11 п— — Я—п— [(пе-1)2 + 2пае-11(пе 11 — к-1)^
+ Я—п—)2 — Я—11п2+п 11 + Я1—11(п2—11-)2 + Я—п—п-1]. (13)
Найдем разностную аппроксимацию для 0е,3/2 , с этой целью сначала проинтегрируем выражение для потока 0е = —Бе дпе/дх — ¡лепеЕ по отрезку [х1, х2]:
Х2 Х2
J Ое ¿X = ! ^ — Бе -X — РепеЕ^ ¿X.
Левую часть равенства (14) аппроксимируем квадратурной формулой цен-
Х 2
тральных прямоугольников У Ое ¿х х Оз/2„, правую разобьем на два интеграла
Х1
+ . Интеграл также аппроксимируем квадратурной формулой центральных прямоугольников, заменив производную разностным отношением
* =,в де) ¿х „—„В
(Х1 +к/2,г3)
1
— Бе(х1 + Н/2,ге)(пе(х2,Ье) — пе(х1 . (15)
Интеграл 3, в свою очередь, разобьем на два интеграла, выделив положительную и отрицательную части Е , получим
Х 2 Х 2
3 = !( — Е) + Мепе ¿х + 1 (—Е)-Мепе ¿х.
Заменим первый из интегралов квадратурной формулой левых прямоугольников, а второй - правых. В результате имеем
35 X Н(—Е) + '3-1 Ме(х1,Ьв )пе(х1,Ьа) + „(—Е)-'"-1 Ме(х2,Ь з)пе(х2,Ь з) .
Далее, согласно методу направленных разностей можно аппроксимировать поток Ое в точке х = 0 одним и тем же (независимо от знака Е) выражением
Се,1/2 = —пе\/ 8кТе,1/(пШе )/4 + Е-^(М+,1'п+,1 + М2++,1).
Действительно, если поле направлено от электрода (то есть Е > 0), то Е1- = 0, а значит,
Се,1/2 = —пе\/ 8кТе,1/(пШе)/4. Поэтому граничное условие (2) при х = 0 заменяется разностным
пв _п3-1
е'1 т е'1 „ — ВЗ-2(п%2 — <д)/„ + ( — Е1/2 — Е3/2)+'3-1ме-1пе,1/2+ + п"еЛ]/ 8кТе,1/(пШе )/4 + ( — Е5/2 — Е3,,2)" м'е-? <,1/2 —
— ( — Е1/2 — Е3/2)-'"-1м1-1п1А/2 = „{я^пе^МГ1 + Щ^п^Ч + щ^п-пе! — я^пеуп+1 + щ-т-)2 — я1-!п2+Хпел + я-п-^ч
+ Н-п^п-' 1 — Н"-11п+ '1((п"е-11)2 + 2пе-11(п"е, 1 — п-))
— 1(Е1/2 + Е3/2)-'3-1 М+-1п+1/2 — 7 (Е1/2 + Е3/2)-' "-1м2+дп2+'1/2. (16)
На нагруженном электроде ( х = Ь ) разностные аналоги граничных условий (3) выводятся интегрированием уравнения (1) по отрезку [Ь — 3Н/2, Ь — Н/2]; так же,
как и в случае граничного условия при х = 0, имеем Н + Ве-М-3/2(. пе, М-1 — пе, М-2)/Н —
"е' М-1 пе' М-1 , , 1
' ' и 1 п® 1 \пе
( — ЕМ-1/2 — ЕМ-3/2) ' М"е- М-1пе, М-1 + К, М-1
2 +4
- ' 3-1 3-1
^8кТе' М-1/(пШе) —
( — ЕМ-5/2 — ЕМ-3/2) ' м5,М-2п&е, М-2 2 '
( — ЕМ-1/2 — ЕМ-3/2) + ' м3е,М-1<, М-1
2
Н"-1 АТ3-1 + Я3-1 (п&-1 )2 + Н8-1 п"-1
я1'М-1' е'М-1ММ-1 "Г я2'М-Ппт,М-1! "Г я3,М-1пт,М-1пе,М-1
— Н4-М-1пе' М-1п+' М-1 — Н5,М-1п+' М-1 [(п1-М-1)2 + 2пе,М-1(пв' М-1 — п1-М-1)] +
+ Н5-1 (п5-1 )2 _ Н5-1 п5 п5 + Н5-1 (п5-1 )2 +
"Г" я10'М-Ппт' М-1! я13'М-1' 2 +,М-1пе,Мя16,М-Пп2*,М-1> '
+ Н5-1 п5-1 п5-1 + 2п5-1 (п5 _ п5-1 ) +
"Г" Н17'М-1пт' М-1п2*' М,М-Ппе,М-1 пе,М-1>\ '
(ЕМ-1/2 + ЕМ-3/2) + ' М+М-1п+ М-1 + 7-2-1-+
(ЕМ-1/2 + ЕМ-3/2)+ ' М2 +' М-1 п2 +' М-1 ,.„4
+ 7-2-. ()
3. Вычисление плотности потока заряженных частиц
Для вычисления потока электронов запишем уравнение (1) в виде системы
дп дС
+ = Я1 пеМ + Я2 п^п + Н3птпе — Н4пеп+ — Я5'п+п2е +
+ Нюп2т — Я13п2+пе + Н1&п\* + Яиптп2*, 0 < х < Ь, Ь> 0, (18) Се = — Ве дпе/дх — пемеЕ. (19)
Вычисление потока проведем, используя модифицированный метод Шарфет-тера-Гуммеля, предложенный для расчета концентрации электронов в полупроводниках в работе [21]. Основным отличием настоящей работы от [21] является использование значений концентрации электронов пе с верхнего временного слоя, найденного при численном решении неоднородного дифференциального уравнения (19) с помощью разностного метода, предложенного в разд. 2 настоящей работы. Кроме того, в отличие от [21], при построении модификации метода Шарфеттера-Гуммеля мы допускаем непостоянство подвижности электронов ме в разряде.
Построим в соответствии с [21] сеточную аппроксимацию аналитического решения пе линейного неоднородного уравнения (19):
С Н 6
пе(х)= [пе, I--е 1+1/2 I ехр(а£') ехр (—а£),
\ Ве,1+1/2 .) )
0
где а = Ме,1 Н1Е1+1/2/Ве,1+1/2, „I = хь+1 — х,, £ = (х — х,)/„, . Тогда при х = х+ , считая концентрации электронов п5 известной, получаем
С5е,1+1/2 = В5-^1/2К, — ехр (а)п^)/(С0 „,), С0 = (ехр (а) — 1)/а, (20)
Необходимое условие устойчивости разностной схемы имеет вид
hi«2Ds-U/2/№-1\AEe-1)\. (21)
Во избежание неустойчивости схемы по аналогии с [21] введем дополнительные точки сетки xl,r = (xi + xi+i)/2 Т К/2, где hu = ^е • 2Dse-+1/2hi/(pse-1\AEs-1\).
Нетрудно убедиться в том, что за счет выбора е G (0,1) можно добиться выполнения условия (21). Значения концентрации во вспомогательных точках ne,L,R будем аппроксимировать следующим образом:
ne,L,R = Kii + 1) exp (a(xi+i Т К/2 - xi)), где a = h-1 log((nlii+i + 1)/nh).
Поскольку мы используем значения концентрации на верхнем временном слое, вычисленные на основе неявной разностной схемы (12), (16), (17) на сетке с постоянным шагом hi = h, то для вычисления плотности потока заряженных частиц мы используем формулу (20), если h < hv .В противном случае, когда h > hv, то, вводя обозначения Wi = 1Esl-1, W+1/2 = -Ve-l\l/2Esl+:l/2, Wi = Wi+i/2 +
+ W—1/2)/2, мы получим следующее выражение для потока (с учетом направленности поля E):
Ge
Ge,i + 1/2
De-1
De,i+1/2
hv G De-1
De,i+1/2
neL - exp(au )(1 + auf3u) neR , \Wl\ > \Wr\,
exp(-a„)(1 + avdv) 1ne,L - ne,R , \Wl\ < \Wr\,
где
Wlr = Wi+i/2 Т AWv/2, AWv = hv(Wi+i + W-)/h, av = -hvWL/Dei + 1/2, Pv = (Wr - Wl)/(2Wl), Sv = hvWR/Dei+1/2, Я = (Wr - Wl)/(2Wr), G =(al - 2avд)/al - exp (-Sv) (a^v + Q-l + 2a.vд + /at,
G = -(al + 2av Pv)/a3 + exp (a.v)(aIfiv + aI + 2av Pv - a^eSj / a3.
Отдельно рассмотрим случай, когда hv > h и E+1/2 = 0. При одновременном выполнении этих условий возникает неопределенность при вычислении потока согласно (20), так как при этом Go = 0. В этом случае воспользуемся предложенной в [21] интерполяцией для Е по формуле ES-1 = ES-1i(1 + 2/Зц), где п = (xi -- xi+i)/h. В результате имеем Д, = AES- ^^ES^), Sv = ^hE'—i/DS-!^ , а значит, поток будет вычисляться следующим образом:
Ge,i+1/2 = De-+1/A exp(-av )(1+ av Д )-1nSej - nSej) / (hv G* ),
где
G* = (a2 - 2avД)/a? - exp (-Sv) (ДЗд + a2v + 2SvД + 2a^v) /a?, то есть G* вычисляется по той же формуле, что и G, но с другими av и /5v.
hv G
В случае, когда \Wl\ = \Wr\ , то есть в ситуации, когда отсутствует конвекция, имеем Gse¡ 1+1/2 = Dse-l1/2(nae , -nse11+1)/h.
Напомним, что, в отличие от метода работы [21], в котором поток вычисляется по концентрации ne с нижнего временного слоя, а значение концентрации на следующем временном слое находится путем численного решения однородного уравнения (nse+1 — nse г)/т = (Gse 1+1/2 — Gse ¡_i/^)/h, в настоящей работе поток рассчитывается по уже найденной концентрации с верхнего временного слоя, что позволяет избежать ограничений на временной шаг, необходимых для выполнения условия устойчивости Куранта, а это значительно уменьшает время расчетов. Результаты численных экспериментов при моделировании ВЧЕ-разряда с применением описанного в настоящей статье метода, полученные в [28], показали хорошее качественно и количественное совпадение с экспериментальными данными [31].
Заключение
ВЧЕ-разряд обладает рядом особенностей, осложняющих его численное моделирование. В частности, в описываемых уравнениями диффузии-конвекции процессах переноса заряженных частиц в сформировавшейся амбиполярной области, расположенной в центре разряда, преобладает диффузия заряженных частиц, а в приэлектродных областях доминирует конвекция. Дополнительное (в силу нелинейности процессов, протекающих в разряде) усложнение вызвано необходимостью определения значений потоков заряженных частиц в ходе решения задач баланса энергии этих частиц. В настоящей работе предложена устойчивая неявная разностная схема первого порядка точности по пространственной переменной для решения этих задач, построенная с помощью интегро-интерполяционного метода (с применением направленных разностей для аппроксимации конвективного слагаемого) и метода Шарфеттера-Гуммеля. Эта схема является двухшаговой; на первом шаге вычисляются значения концентрации заряженных частиц на верхнем временном слое, а на втором шаге - значения плотности потока заряженных частиц по уже найденной концентрации. Использование концентрации с верхнего временного слоя позволяет избежать ограничений, вызванных необходимостью удовлетворения условию Куранта, на временной шаг, и, следовательно, уменьшить время вычислений. Применение построенной модификации метода Шарфеттера-Гуммеля уменьшает погрешность при вычислении потока.
Благодарности. Работа выполнена при финансовой поддержке РНФ в рамках научного проекта 16-11-10299.
Литература
1. Полак Л.С., Синярев Г.Б., Словецкий Д.И. и др. Химия плазмы / Под ред. Л.С. По-лака, Ю.А. Лебедева. - Новосибирск: Наука, 1991. - 325 с.
2. Бердичевский М.Г., Марусин В.В. Нанесение покрытий, травление и модифицирование полимеров с использованием низкоэнтальпийной неравновесной плазмы: Обзор. -Новосибирск: Ин-т теплофизики, 1993. - 107 с.
3. Gadzhiev M.K., Isakaev E.K., Tyuftyaev A.S., Yusupov D.I. A high-power low-temperature air plasma generator with a divergent channel of the output electrode // Techn. Phys. Lett. - 2016. - V. 42, No 1. - P. 79-81. - doi: 10.1134/S106378501601020X.
4. Бадриев И.Б., Желтухин В.С., Чебакова В.Ю. О решении некоторых нелинейных краевых и начально-краевых задач // Материалы XXII Междунар. симпозиума
«Динамические и технологические проблемы механики конструкций и сплошных сред» им. А.Г. Горшкова. - М.: ООО «ТРП», 2016. - С. 31-33.
5. Каренгин А.Г. Физика и техника низкотемпературной плазмы. - Томск: Изд-во Том. политехн. ун-та, 2008. - 129 с.
6. Lebedev Yu.A., Tatarinov A.V. Electrodynamics of microwaves in a coaxial non-regular waveguide partly filled with plasma // Plasma Sources Sci. Technol. - 2004. - V. 13, No 1. - P. 1-7. - doi: 10.1088/0963-0252/13/1/001.
7. Амиров Р.Х., Исакаев Э.Х., Шавелкина М.Б., Юсупов Д.И., Шаталова Т.Б., Эмиров Р.М. Плазмоструйный синтез углеродных наноструктур и устройство для его реализации // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2014. - Т. 156, кн. 4. - С. 112-119.
8. Лебедев Ю.А., Та,та,ринов А.В., Титов А.Ю., Эпштейн И.Л. Двумерная модель неравновесного сильно неоднородного СВЧ-разряда во внешнем постоянном поле // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2014. - Т. 156, кн. 4. - С. 120-132.
9. Wang H.-B, Sun W.-T., Lia H.-P., Bao Ch.-Y. Discharge characteristics of atmospheric-pressure radio-frequency glow discharges with argon/nitrogen // Appl. Phys. Lett. -2006. - V. 89, No 16. - Art. 161504, P. 161504-1-161504-3. - doi: 10.1063/1.2362631.
10. Лебедев Ю.А. О сравнении неравновесной плазмы самостоятельных электрических разрядов // Высокочастотный разряд в волновых полях: Сб. науч. тр. - Самара, 1990. - С. 25-63.
11. Абдуллин И.Ш., Желтухин В.С., Чебакова В.Ю. Высокочастотный емкостной разряд: моделирование (обзор) // Вестн. Казан. технол. ун-та. - 2014. - Т. 17, № 23. -С. 9-14.
12. Тихонова Н.В., Желтухин В.С., Чебакова В.Ю., Бородаев И.А. Математическая модель высокочастотной плазменной обработки многослойных материалов заготовки верха обуви // Вестн. Казан. технол. ун-та. - 2012. - Т. 15, № 17. - С. 36-39.
13. Goedheer W.J. Lecture notes on radio-frequency discharges, dc potentials, ion and electron energy distributions // Plasma Sources Sci. Technol. - 2000. - V. 9. - P. 507516. - doi: 10.1088/0963-0252/9/4/306.
14. Райзер Ю.П., Шнейдер М.Н., Яценко Н.А. Высокочастотный емкостный разряд: Физика. Техника эксперимента. Приложения. - М.: Изд-во МФТИ, 1995. - 320 с.
15. Райзер Ю.П. Физика газового разряда. - Долгопрудный: Изд. дом «Интеллект», 2009. - 736 с.
16. Chebert P., Braithwaite N. Physics of radio-frequency plasmas. - Cambridge: Cambridge Univ. Press, 2011. - 386 p.
17. Бадриев И.Б., Чебакова В.Ю. Математическое моделирование низкотемпературной ВЧЕ-плазмы в аргоне // Математическое моделирование и краевые задачи: Труды Десятой Всерос. науч. конф. с междунар. участием: в 3 ч. - Самара: СамГТУ, 2016. -Ч. 2. - С. 17-21.
18. Badriev I.B., Chebakova V.Y., Zheltukhin V.S. Capacitive coupled RF discharge: modelling at the local statement of the problem // J. Phys.: Conf. Ser. - 2017. - V. 789, No 1. - Art. 012004, P. 1-4. - doi: 10.1088/1742-6596/789/1/012004.
19. Бадриев И.Б., Желтухин В.С., Чебакова В.Ю. Численное исследование низкотемпературной ВЧЕ-плазмы в аргоне при повышенных давлениях // Материалы XI Междунар. конф. по неравновесным процессам в соплах и струях (NPNJ'2016). - М.: Моск. авиац. ин-т, 2016. - С. 477-479.
20. Scharfetter D.L., Gummel H.K. Large-signal analysis of a silicon Read diode oscillator// IEEE Trans. Electron Devices. - 1969. - V. 16, No 1. - P. 64-77. - doi: 10.1109/T-ED.1969.16566.
21. Kulikovsky A.A. A more accurate Scharfetter-Gummel algorithm of electron transport for semiconductor and gas discharge simulation //J. Comput. Phys. - 1995. - V. 119, No 1. - P. 149-155. - doi: 10.1006/jcph.1995.1123.
22. Дегтярев Л.М., Фаворский А.П. Потоковый вариант метода прогонки для разностных задач с сильно меняющимися коэффициентами // Журн. вычисл. матем. и ма-тем. физики. - 1969. - Т. 9, № 1. - С. 211-218.
23. Гольштейн Е.Г., Третьяков Н.В. Модифицированные функции Лагранжа. - М.: Наука, 1989. - 400 с.
24. Бадриев И.Б., Задворнов О.А. Итерационные методы решения вариационных неравенств второго рода с обратно сильно монотонными операторами // Изв. вузов. Матем. - 2003. - № 1. - С. 20-28.
25. Бадриев И.Б., Задворнов О.А. О сходимости итерационного метода двойственного типа решения смешанных вариационных неравенств // Дифференц. уравнения. -2006. - Т. 42, № 8. - С. 1115-1122.
26. Glowinski R., Le Tallec P. Augmented Lagrangian and operator-splitting methods in nonlinear mechanics. - Philadelphia, PA: SIAM , 1989. - x, 292 p.
27. Чебакова В.Ю. Численное моделирование высокочастотного емкостного разряда // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2015. - Т. 157, кн. 2. - С. 126-140.
28. Чебакова В.Ю. Моделирование высокочастотного емкостного разряда при атмосферном давлении в аргоне // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2016. -Т. 158, кн. 3. - С. 404-423.
29. Самарский А.А., Гулин А.В. Численные методы. - М: Наука, 1989. - 432 с.
30. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. - М.: Эдиториал УРСС, 1999. - 248 с.
31. Лисовский В.А. Особенности а - y-перехода в ВЧ разряде низкого давления в аргоне // Журн. техн. физики. - 1998. - Т. 68, Вып. 5. - С. 52-60.
Поступила в редакцию 15.01.17
^Келтухин Виктор Семенович, доктор физико-математических наук, профессор кафедры плазмохимических и нанотехнологий высокомолекулярных материалов Казанский национальный исследовательский технологический университет
ул. К. Маркса, д. 68, г. Казань, 420015, Россия E-mail: vzheltukhin@gmail.com
Фадеева Мария Сергеевна, студент Института вычислительной математики и информационных технологий
Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: manysha-98@mail.ru
Чебакова Виолетта Юрьевна, кандидат физико-математических наук, доцент кафедры анализа данных и исследования операций
Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: vchebakova@mail.ru
МОrЦH®HКАЦHa METO^A mAPOETTEPA-rYMME^a
455
ISSN 2541-7746 (Print)
ISSN 2500-2198 (Online)
UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI
(Proceedings of Kazan University. Physics and Mathematics Series)
2017, vol. 159, no. 4, pp. 444-457
Modification of the Scharfetter—Gummel Method for Calculating the Flux of Charged Particles for Simulation of a Radio-Frequency Capacitive Coupled Discharge
V.S. Zheltukhina* , M.S. Fadeevab** , V.Ju. Chebakovab***
aKazan National Research Technological University, Kazan, 420015 Russia bKazan Federal University, Kazan, 420008 Russia E-mail: *vzheltukhin@gmail.com,, **manysha-98@mail.ru, ***vchebakova@mail.ru
Received January 15, 2017 Abstract
A numerical method of solving of nonstationary convection-diffusion equation for charged particle densities in capacitive coupled RF discharge is presented. The method allows us to calculate the particle concentration together with the flux density simultaneously at very rapid changes of particle density and electrical field in the electrode sheath. The method is a modification of the Scharfetter-Gummel algorithm. Firstly, implicit difference approximation of the equation constructed by the integro-interpolation method is used for calculating the charged particles density. Then the flux density of charged particles is calculated. Our modification allows to carry out the calculation at milder restriction on the time step than the Courant condition. In addition, diffusion coefficient may be variable unlike with the original Scharfetter-Gummel algorithm.
Keywords: mathematical modeling, radio-frequency capacitive coupled discharge, modified Scharfetter-Gummel method, convection-diffusion equation, flux density of charged particles
Acknowledgments. The study was supported by the Russian Science Foundation (project no. 16-11-10299).
References
1. Polak L.S., Sinyarev G.B., Slovetsky D.I. et al. Plasma Chemistry. Polak L.S., Lebe-dev Yu.A. (Eds.). Novosibirsk, Nauka, 1991. 325 p. (In Russian)
2. Berdichevskii M.G., Marusin V.V. Coating, Etching, and Modifying Polymers Using Low-Enthalpy Non-Equilibrium Plasma: Review. Novosibirsk, Inst. Teplofiz., 1993. 107 p. (In Russian)
3. Gadzhiev M.K., Isakaev E.K., Tyuftyaev A.S., Yusupov D.I. A high-power low-temperature air plasma generator with a divergent channel of the output electrode. Tech. Phys. Lett., 2016, vol. 42, no. 1, pp. 79-81. doi: 10.1134/S106378501601020X.
4. Badriev I.B., Zheltukhin V.S., Chebakova V.Yu. On solving of some nonlinear boundary and initial boundary value problems. Materialy XXII Mezhdunar. simpoziuma "Dinamicheskie i tekhnologicheskie problemy mekhaniki konstruktsii i sploshnykh sred" im. A.G. Gorshkova [Proc. XXII Int. Symp. "Dynamic and Technological Problems
of Mechanics of Constructions and Continuous Media" Dedicated to A.G. Gorshkov]. Moscow, TRP, 2016, pp. 31-33. (In Russian)
5. Karengin A.G. Physics and Technology of Low-Temperature Plasma. Tomsk, Izd. Tomsk. Politekh. Univ., 2008. 129 p. (In Russian)
6. Lebedev Yu.A., Tatarinov A.V. Electrodynamics of microwaves in a coaxial non-regular waveguide partly filled with plasma. Plasma Sources Sci. Technol., 2004, vol. 13, no. 1, pp. 1-7. doi: 10.1088/0963-0252/13/1/001.
7. Amirov R.Kh., Isakaev E.Kh., Shavelkina M.B., Yusupov D.I., Shatalova T.B., Emirov R.M. Plasma-jet synthesis of carbon nanostructures and a device for its implementation. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2014, vol. 156, no. 4, pp. 112-119. (In Russian)
8. Lebedev Yu.A., Tatarinov A.V., Titov A.Yu., Epstein I.L. Two-dimensional model of a non-equilibrium strongly non-uniform microwave discharge in a DC external field. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2014, vol. 156, no. 4, pp. 120-132. (In Russian)
9. Wang H.-B., Sun W.-T., Lia H.-P., Bao Ch.-Y. Discharge characteristics of atmospheric-pressure radio-frequency glow discharges with argon/nitrogen. Appl. Phys. Lett., 2006, vol. 89, no. 16, art. 161504, pp. 161504-1-161504-3. doi: 10.1063/1.2362631.
10. Lebedev Yu.A. Comparison of the Non-Equilibrium Plasma of Self-Maintained Electrical Discharges. Radio-Frequency Discharge in Wave Fields. Samara, 1990, pp. 25-62. (In Russian)
11. Abdullin I.Sh., Zheltukhin V.S., Chebakova V.Yu. Radio-frequency capacitive coupled discharge: Modeling (review). Vestn. Kazan. Tekhnol. Univ., 2014, vol. 17, no. 23, pp. 914. (In Russian)
12. Tikhonova N.V., Zheltukhin V.S., Chebakova V.Yu., Borodaev I.A. Mathematical model of radio frequency plasma processing of multilayer materials of shoe uppers. Vestn. Kazan. Tekhnol. Univ., 2012, vol.15, no.17, pp. 36-39. (In Russian)
13. Goedheer W.J. Lecture notes on radio-frequency discharges, dc potentials, ion and electron energy distributions. Plasma Sources Sci. Technol., 2000, vol. 9, pp. 507-516. doi: 10.1088/0963-0252/9/4/306.
14. Raizer Y.P., Shneider M.N., Yatsenko N.A. Radio-Frequency Capacitive Discharges. Boca Raton, Fla., CRC Press, 1995. 304 p.
15. Raizer Yu.P. Gas Discharge Physics. Dolgoprudnyi, Intellekt, 2009. 736 p. (In Russian)
16. Chebert P., Braithwaite N. Physics of Radio-Frequency Plasmas. Cambridge, Cambridge Univ. Press, 2011. 386 p.
17. Badriev I.B., Chebakova V.Yu. Mathematical simulation of low-temperature RF plasma in argon. Matematicheskoe modelirovanie i kraevye zadachi: Trudy Desyatoi Vseros. nauch. konf. s mezhdunar. uchastiem [Mathematical Modeling and Boundary Value Problems: Proc. 10th All-Russ. Sci. Conf. with Int. Participation]. Vol. 2. Samara, 2016, pp. 17-21. (In Russian)
18. Badriev I.B., Chebakova V.Y., Zheltukhin V.S. Capacitive coupled RF discharge: Modelling at the local statement of the problem. J. Phys.: Conf. Ser., 2017, vol. 789, no. 1, art. 012004, pp. 1-4. doi: 10.1088/1742-6596/789/1/012004.
19. Badriev I.B., Zheltukhin V.S., Chebakova V.Yu. Numerical study of low-temperature RF plasma in argon at elevated pressures. Materialy XI Mezhdunar. konf. po neravnovesnym, protsessam v soplakh i struyakh (NPNJ'2016) [Proc. XI Int. Conf. on Non-Equilibrium Processes in Nozzles and Jets (NPNJ'2016)]. Moscow, 2016, pp. 477-479. (In Russian)
20. Scharfetter D.L., Gummel H.K. Large-signal analysis of a silicon Read diode oscillator. IEEE Trans. Electron Devices, 1969, vol. 16, no. 1, pp. 64-77. doi: 10.1109/T-ED.1969.16566.
21. Kulikovsky A.A. A more accurate Scharfetter-Gummel algorithm of electron transport for semiconductor and gas discharge simulation. J. Comput. Phys., 1995, vol. 119, no. 1, pp. 149-155. doi: 10.1006/jcph.1995.1123.
22. Degtyarev L.M., Favorskii A.P. Flow variant of the sweep method for difference problems with strongly varying coefficients. USSR Comput. Math. Math. Phys., 1969, vol. 9, no. 1, pp. 285-294. doi: 10.1016/0041-5553(69)90021-4.
23. Golshtein E.G., Tretyakov N.V. Modified Lagrangians. Moscow, Nauka, 1989. 400 p. (In Russian)
24. Badriev I.B., Zadvornov O.A. Iterative methods for solving variational inequalities of the second kind with inversely strongly monotone operators. Russ. Math., 2003, vol. 47, no 1, pp. 18-26.
25. Badriev I.B., Zadvornov O.A. On the convergence of the dual-type iterative method for mixed variational inequalities. Differ. Equations, 2006, vol. 42, no. 8, pp. 1180-1188. doi: 10.1134/S001226610608012X.
26. Glowinski R., Le Tallec P. Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics. Philadelphia, PA, SIAM, 1989. x, 292 p.
27. Chebakova V.Yu. Numerical simulation of the high-frequency capacitive discharge. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2015, vol. 157, no. 2, pp. 126-140. (In Russian)
28. Chebakova V.Yu. Modeling of radio-frequency capacitive discharge under atmospheric pressure in argon. Lobachevskii J. Math., 2017, vol. 38, no. 6, pp. 1165-1178. doi: 10.1134/S1995080217060154.
29. Samarskii A.A., Gulin A.V. Numerical Methods. Moscow, Fizmatlit, 1989. 432 p. (In Russian)
30. Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Convection-Diffusion Problems. Moscow, Editorial URSS, 1999. 248 p. (In Russian)
31. Lisovskii V.A. Features of the a - 7 transition in a low-pressure rf argon discharge. Tech. Phys., 1998, vol. 43, no. 5, pp. 526-534. doi: 10.1134/1.1259033.
Для цитирования: Желтухин В.С., Фадеева М.С., Чебакова В.Ю. Модификация / метода Шарфеттера-Гуммеля для нахождения потока заряженных частиц при мо-\ делировании высокочастотного емкостного разряда // Учен. зап. Казан. ун-та. Сер.
Физ.-матем. науки. - 2017. - Т. 159, кн. 4. - С. 444-457.
For citation: Zheltukhin V.S., Fadeeva M.S., Chebakova V.Ju. Modification of the / Scharfetter-Gummel method for calculating the flux of charged particles for simulation of \ a radio-frequency capacitive coupled discharge. Uchenye Zapiski Kazanskogo Universiteta.
Seriya Fiziko-Matematicheskie Nauki, 2017, vol. 159, no. 4, pp. 444-457. (In Russian)
