Об одной численной схеме экспоненциальной подгонки для решения уравнений высокочастотного разряда в гидродинамическом приближении Текст научной статьи по специальности «Физика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2014 Математика и механика № 3(29)

УДК 533:519.6

А.Г. Горобчук

ОБ ОДНОЙ ЧИСЛЕННОЙ СХЕМЕ ЭКСПОНЕНЦИАЛЬНОЙ ПОДГОНКИ ДЛЯ РЕШЕНИЯ УРАВНЕНИЙ ВЫСОКОЧАСТОТНОГО РАЗРЯДА В ГИДРОДИНАМИЧЕСКОМ ПРИБЛИЖЕНИИ1

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

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

Исследование физики ВЧ-разряда является сложной многопараметрической задачей, которая носит преимущественно исследовательский характер. В диапазоне рабочего давления 0,1—3,0 торр относительно полное представление о распределении концентраций электронов и ионов, а также температуры электронов дает моделирование ВЧ-разряда в гидродинамическом приближении [1, 2]. Умеренная вычислительная трудоемкость гидродинамического подхода позволяет применить его для адекватного исследования процессов, реализуемых при плазмохимическом травлении полупроводниковых материалов.

На стадии горения ВЧ-разряда, характеризующейся образованием тонких при-электродных слоев и квазинейтрального столба, из-за сильных градиентов параметров плазмы и большой проводимости проведение расчетов с применением явных численных схем для решения уравнений непрерывности электронов и ионов становится крайне неэкономичным [3]. Из-за ограничений, накладываемых условием Куранта для электронов Ке = Ует / к < 1 (где Уе - дрейфовая скорость электронов, т, к - шаги по времени и по пространству), шаг по времени оказывается много меньше, чем характерное время развития разряда, что может потребовать более 104 - 105 шагов по времени. Здесь эффективно применение неявных конечно-разностных схем, позволяющих использовать шаг по времени в 10 - 102 раз больший, чем дает условие Куранта. Однако в численных схемах, в которых распределение потенциала находится непосредственно из уравнения Пуассона, возникает дополнительное ограничение вида т< 1 / 4пст, где ст - проводимость плазмы [3]. Данное ограничение дает значение на порядок больше, чем условие Куранта. Кроме того, основные уравнения обычно решаются с использованием дополнительных итераций самосогласования. Само решение системы нелинейных разностных уравнений также представляет большую трудность. В настоящее время при моделировании газовой плазмы напряженность электрического поля находится из уравнения сохранения электрического тока. Для решения данного уравнения применяются консервативные по заряду разностные схемы, в которых шаг по времени ограничивается периодом лэнгмюровских колебаний ионов. При использовании чисто неявных схем в уравнении тока решение разностных уравне-

1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект №11-01-00064), президентской программы для государственной поддержки ведущей научной школы РФ (проект НШ-6293.2012.9).

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

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

Физико-математическая модель

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

Введем переменные:

П ~ М-1 ^ В1 ,

р 1 _ —, м _—, 1^1 _ —,1 _е, р;

Пе0 М10 В10

ф п Те Г 2 Ь

Ф = —, ^е =Т±-, т = Я %= -, С = Ь , А _ -.

ф0 те0 Г Ь Г

Здесь индексы е, р обозначают электроны и ионы соответственно; щ - плотность частиц; ц1, Б1 - подвижность и коэффициент диффузии; Те - электронная температура; ф - потенциал; / - частота активации; t - время. Плотность, подвижность, коэффициент диффузии, потенциал и температура обезразмеривались на среднюю плотность электронов пе0, характерные значения ц10, Б10, потенциал ВЧ-элек-

трода Ф0 и температуру Те0, соответствующую тепловой энергии электронов. Координаты нормированы на радиус цилиндра г и межэлектродное расстояние Ь соответственно.

Концентрации плазменных компонент находятся из уравнений непрерывности для электронов и ионов:

др1 + . + Л} А_«с Т_ Ь п 1, 1 _ е

Р1 Т-дТ- + А дт

■ + -

д£

+ —±_Р,8,, Т_^~, р1 Н , (1)

дС Ве0 1 Щ/ Вр, 1 _ р.

Потоки частиц определяются формулами

„ _ дф дВ 1 р1 _ дф дВ 1 р1

]) _-*•)Рб/М-,р,-, ]) _-в,Ре,ц,р,--------------------------

% 1 111 д| д% ^ 1 д^ дй,

Ре, =

Н-гоф о О,0

51 ={-11,

I = е, I = р.

Нижние индексы | и ^ относятся к направлениям радиальной и аксиальной

координат. Число Пекле характеризует соотношение процессов переноса и диффузии. Источниковый член включает ионизацию исходных молекул и процесс прилипания электронов к ионам:

=

{Ре (Ба Г - га X I Ре Ба Г,

I = е, 1 = p,

Г = ехр

Ве,

Ба =-

I2 пк,0

Ом

-ехр

1 -±

Е

кТеГ

А

кГ„,

уе0 е0 V ГУ1е0,

где г - скорость ионизации, га - скорость прилипания, Е - энергия ионизации, к - постоянная Больцмана, va - частота прилипания, п - объемная плотность газа, к0 - коэффициент скорости ионизации. Число Дамкелера пропорционально отношению интенсивности ионизации к диффузии.

Распределение электрического потенциала определяется из уравнения Пуассона:

,2 1 д дф%

’д|.

А_

+ |_ = Г(Ре-Рр ), Г =

єФ0

I д| V ^ / (Л, ь^0

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

Электронная температура находится из уравнения баланса электронной энергии:

д , „ ч ,21К + д

3

- Тд-(ре 9е)+ А2

5 дт

35

(3)

где компоненты потока электронной энергии имеют вид _ ]е^е, _ Уе?^е.

Параметр х характеризует отношение энергии электрического поля и тепловой энергии электронов. Другой параметр © отвечает за потери энергии в реакциях ударной ионизации. Их выражения имеют вид

еФ0

Х ='

К,

© = ■

к

к.

ке0 = _ kTe0,

е0 е0

где ке0 - тепловая энергия электронов, к - потери энергии при ионизации.

На границах расчетной области для уравнений (1) - (3) ставятся следующие краевые условия:

- на нижнем электроде (0 <5<5о, С = 0):

дО р р р

5С

= 0, ф = 0, 9е =9е

- на верхнем электроде (0 <5<5о, С = 1):

= 0, ф = 5Іп(2лт), 9е = 9е

Ре = -и

дО р р р

5С

- на оси симметрии (5 = 0, 0 < С < 1):

= 0, ^ = 0, ^ф = 0, ^ = 0;

55 55 55 55

- на боковой поверхности (5 = 1, 0 < С < 1):

Ре = 0, 90 рР р = 0, ^5 = 0, 9е =9е5;

г е ' ' е е5 '

55 55

- на перфорированной границе откачки (5о < 5 < 1, С = 0):

дВ р р р Зф

Ре = 0, = 0, ^ = 0, 9е = ^;

5^ 5^

- на верхнем электроде (50 < 5 < 1, С = 1):

дВ р р р дф>

Ре = 0, ---р-^ = 0, ^ = 0, 9е =9е5;

е 5С 5С е

где 5о = го / г - радиус электродов, у - коэффициент вторичной электронной

эмиссии, 9е5 - электронная температура на электродах.

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

Конечно-разностная схема

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

Основным методологическим подходом при построении безусловно монотонных разностных схем является их регуляризация, частным случаем которой является метод экспоненциальной подгонки или интегральных тождеств, использованный в данной работе [4]. Рассмотрим суть этого метода на примере уравнений непрерывности (1), записанных для полных потоков частиц, включающих одновременно конвективную и диффузионную части. В рассматриваемой стандартной области введем прямоугольную сетку с узлами ^/), / = 1,...,N; / = 1,...,М. Вся

область разбита на некоторое количество непересекающихся контрольных объемов. Проинтегрируем дифференциальное уравнение (1) по выделенному контрольному объему. Первоначально преобразуем выражение для компоненты потока / к следующему виду:

Л

- + О і Рі - -

т. Ь дф ■с ■с

где ыг = Ре, —-------, и проинтегрируем его в пределах от ^ до §/+1:

Б1

1 і Рі

Л

і рі

--(иі ?)1!Г •

Предполагается, что потоки частиц, скорости дрейфа и коэффициенты переноса

постоянны между узлами разностной сетки: Б,р, I = £), р, , Б,р, I = £), р, ,

/ / ^/+1 '+1 г'+1

Лі,

- Л , и I?'+1 - и •^■+1/2 і іі+

• Таким образом, приходим к равенству

ґ

1п

-+В і Рі

іі +1 г іі +1

ґ

- 1п

4.

-+0 і Рі,.

+ 1/2 Ь/ +1/2 •

V / +1/2 / V ,/ +1/2

где Д^+ = 5/+1 -§/. После преобразования получаем выражение для потока

частиц

Р,+1 А ехр(м, 1/2 Д5/+1/2) -р, £>,

_ /+1 /+1 /+1/2 ^+1/2 / /

Л

?/+1/2

і

ехр(и

.+1/2 5/+1/

Данное равенство справедливо для любого / = 1,...,М , поэтому окончательное выражение имеет вид

Д5 ) — Р,

Лі,

?/+1/ 2, у

- -и

Рі ° ехр(иі Д? )-Рі Ві

/+ 1,у /+1,у /+1/2, у Ъ/+1/2 /, у / ,у

+1/2,л ехр(иі Д? ) -1

і/ +1/2, у '+1/2

Аналогичным образом вычисляются остальные потоки / , / , / .

-1/2, _/ , _/ + 1/2 ^/, _/-1/2

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

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

к+1 к к+1 'к+1 'к+1 . -к+1 -к+1 -к+1

Р; - Р; ЛЕ - А ЛЕ + ЛЕ ]гс - А

р т 1,J_____+ А +1/2’J_^~1/2,, +_Ь/+1/2^_____./-1/2,, _____5/,J-1/2 _ р £к

Ч Д. 2. ДС

где потоки частиц в радиальном и аксиальном направлениях вычислялись по

формулам

2к+1 ^к+1 т\к ^..,.,.4/ „к+1 \ ^к+1 7^к

■к+1 _ ^/ +1/2,, ' */ +1,, */ + 1,,

Д '

Л+1 7^ к Лжт/ „к+1 \ +1 Л к

Р1 °1 ехР (2 -Р1 °1

> /+1,, /+1, ■_______\ / +1/2, , / /, , /,

./+1/2 ехр (2^ )- 1

к+1 к+1 к ( к+1 \ к+1 7^ к

2 Р1 А ехР (2 )-Р1 £>1

■к+1 _ _ /,■+1/2 /,j+1 /,■+1 V /,,+ 1/2 / /,, /,,

, *

./+1/2, , д.

+1/2

' - 2,,

С,..1/2 ехр(■2 )- 1

Здесь индекс к относится к моменту времени тк, при этом временной шаг составляет величину Дт _ тк+1 - тк . Символы 2<к+1 и 2^+1 даются формулами

к /+1/2,, /,,+1/2

дк дк

2к _-Лре ^+1/2,, к -шк) 2к _-лРе ^+1/2 к -шк)

21/ +1/2,, _ ^ +1,, ф,,), 2/,, +1/2 _ ^ ) 1ф,,+1 Ф.',, ) ■

'/ +1/2,, '/,, + 1/2

В направлении . узел /+1/2 располагается посередине между узлами / и /+1 разносной сетки, которые разделены отрезком Д.. _ .+1 -.. Аналогичное распо-

ложение узлов используется в направлении ^ , где узел ,+1/2 располагается посередине между узлами , и ,+1, разделенных отрезком Д^, 1/2 _С,+1 - С, ■ Другие

обозначения имеют вид Д.. _(.+1 -.-1)/2, Д^, _(^,+1 -С,-1)/2. Выражения

для потоков ,+1 , ,,к+1 в узлах сетки (/-1/2,,) и (/, ,-1/2) имеют аналогич-

'.-1/2,, 'С/,,-1/2

ное представление. Значения потенциала фк+1,, фк,, фк,+1, коэффициентов пек к к к к реноса Д' , Д' , ) , ) и источниковый член О' в численной

г г '/+1/2,, Г '/,,+ 1/2 '/+1/2,, '/,,+1/2 '

схеме брались с предыдущего временного слоя тк.

Для нахождения электрического потенциала решается уравнение Пуассона с помощью неявной конечно-разностной схемы типа «крест»:

( к к к к \ ^ ( к к к к \

А

'Т - й&- - &- 1 ф/,,+1 -ф/,, ф/,, -ф /,,-1

д 1ГК 1+ 1 д 1ГК 1дс, д + 1 д

=ГК Ч., )■

Аналогичным образом интегрируются выражения для компонент потока электронной энергии. Значения электронной температуры в узлах сетки находятся из уравнения

3 Зк+у+1 -9к „к 4+1 - ,к+1 ,*+1 + 4+1

3 ^е,,,Ре,,, е,,,Ре,-,, + А 2 +1/2,, ^./-1/ 2,, + А 2 . +1/2,, ./-1/2,, ,

5 ДТк Д., 2.-

■£+1 -£+1

JSr , JSr

+ ri,J+1/2 ri>J-'

—+x a2 (j£+1 e£+1 + j£+1 e£+1 ) +

2 \ % +1/2, J i +1/2, J Je^i-l/2,j i-1/2, j)

+x( j£+1 eZ+1 + j£+1 eZ+1 )+©s£ = 0,

2v eZi,j+1/2 Zi,j+1/2 %■,j-1/2 ^j-i/2/ ^j

где компоненты потока электронной энергии в радиальном и аксиальном направлениях вычисляются по формулам

£ z£+l s!+1 pe+1 d£ exp (z£+1 )-s!+1 р!+14"

•£+1 i +1/2,J i +1,J i+1,J i+1,J V i +1/2,J f i,J i,J i,J

^ +1/2, J Д

" +1/2

.i+1/2 exp (z£+1 )-1

V i +1/2, J /

z!+1 S+1 p!+1 d£ exp (zke +1 )-»e+1p!+1 Dke

i,J +1 i,J+1 \ i,J+1/2 / i,J i,J i

j£+1 = i, J+1/2 i, J+1 i, J +1 i, J +1 V i, J+ 1/2 / i, J i, J i, J

JSz _ A

Zi, J +1/2 Д.

/2

4jexp(J )- 1

Компоненты электрического поля определяются как

т?£+1 /,_£+1 ,_£+1 \ / а 7т,£+1 (,Jc+1 £+1 \ / а

Es i+1/2,j =-((P<+W-Ъ.1 )7Д.,+1/2, EC„^+1/2 = -(cJ -^J )'ДсJ+1/2-

Выражения для компонент потока электронной энергии jS+1 , jS+1 и элек-

%i-1/2, J zi, J-1/2

трического поля E£+1 , EZ+1 имеют аналогичное представление. Вследствие

^i-1/2, J i, J-1/2

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

Учитывая наличие больших градиентов электрического поля, для устойчивого расчета потоков электронов и ионов на границах области (на электродах и боковой стенке камеры) использовалась формула линейной интерполяции, позволяющая выразить поток на границе через значения потока внутри области, рассчитанного по экспоненциальной схеме, например, на верхнем электроде:

a j£+1 _ j£+1

Дг

:£ + 1 = ■£ + 1 , z J +1/2 Zi, J +1/2 Zi, J-1/2

J Sz 1 / О Л A -

Zi,j+1 Zi,j+1/2 2 Дг

Z J

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

Результаты и обсуждение

В расчетах моделировался ВЧ-разряд в плазмохимическом реакторе травления радиальной схемы при рабочем давлении p = 0,5 торр [5]. Радиус реакционной камеры выбирался равным rt = 30 см, а межэлектродное расстояние L = 3,5 см.

ВЧ-разряд инициировался на рабочей частоте f = 13,56 МГц при напряжении на электродах Ф0 = 110 В. Другие параметры выбирались следующими: плотность газа n = 3-1016 см-3, средняя плотность электронов ne0 = 6-109 см-3, электронная температура на электродах и боковой стенке Tes = 0,5 эВ. Другие размерные величины, входящие в уравнения, брались из работ [1, 2, 5, 6]. Начальные распределения концентраций электронов и ионов выбирались в виде аналитического решения, полученного в диффузионно-дрейфовом приближении для ВЧ-разряда в цилиндрической области [7]: ре0 =р 0 = 3,63796 J0 (2,405^) sin(nZ). Начальное

распределение электронной температуры полагалось постоянным (рис. 1).

Рис. 1. Изолинии электронной плотности ne -10 10 см 3 в плазмохимическом реакторе

Моделирование ВЧ-разряда с гармоническим током разряда требует решения нестационарных уравнений переноса зарядов и баланса энергии электронов. В расчетах рассматривалось до 104 периодов ВЧ-разряда, когда колебания электронной плотности и температуры выходили на периодический режим. Из-за высокой нелинейности исходных уравнений, включающих решение уравнения Пуассона, полученные неявные численные схемы являются условно устойчивыми. В этой связи шаг по времени выбирался в диапазоне Дт » 5 • 10-6 -10-4 , при этом

размер сетки составлял 40x20 узлов по пространственным переменным.

Тестирование алгоритма проводилось на основе сравнения численного решения с аналитическим, полученным в диффузионно-дрейфовом приближении для ВЧ-разряда между двух поверхностей [7]: ре0 = 1,57079Бт (л^). При этом на

верхнем электроде (0 ^5о, С = 1) вместо краевого условия ре = -/. бралось

условие ре = 0 , а на выходе из плазменной области (£ = 1о,0<С< 1) ставились мягкие условия. Найдено, что расчетная электронная плотность однородна вдоль радиальной координаты, а ее профиль на оси симметрии практически совпадает с синусоидальной зависимостью от высоты реактора. Максимальное локальное отклонение не превосходит величины 0,05, что находится в пределах ошибки аппроксимации.

На рисунке представлено распределение электронной плотности в плазмохимическом реакторе травления радиальной схемы. Электронная плотность максимальна в центре реактора и монотонно спадает до нуля с приближением к электродам, а также и в радиальном направлении к боковой стенке реактора. Максимальная плотность электронов (в сечении г = 3,25 см) монотонно уменьшается вдоль радиуса на 0,3, 2,3, 17,1 и 72 % при г = 8,14, 14,57, 21,00 и 27,43 см соответственно. Если не рассматривать крайнее значение г = 27,43 см, то снижение расчетного распределения электронной плотности по радиусу реактора составляет

~17 %, что должно существенно повлиять на однородность обработки кремниевых подложек. Отношение безразмерных времен для электронов и ионов составляет величину 104, что соответствует более быстрому процессу диффузии электронов по сравнению с ионами. В расчетах это приводит к тому, что ионы остаются неподвижными, а электроны откликаются на высокочастотное электрическое поле. Сравнение расчетной электронной плотности в различные моменты периода колебаний напряжения показали, что концентрация электронов изменяется только в близи электродов. Распределение ионной плотности повторяет структуру распределения электронной плотности и не реагируют на высокочастотное поле. Среднее по периоду значение электронной температуры приблизительно равно 5,9 эВ. Исключение составляют узкие приэлектродные слои, где электронная температура резко меняется.

Заключение

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

ЛИТЕРАТУРА

1. Graves D.B., Jensen K.F. A continuum model of DC and RF discharges // IEEE Transactions on Plasma Science. 1986. V. PS-14. P. 78-91.

2. Lymberopoulos D.P., Economou D.J. Fluid simulation of glow discharges: effect of metastable atoms in argon // J. Appl. Phys. 1993. V. 73. P. 3668-3679.

3. Гадияк Г.В., Швейгерт В.А., Ууэмаа О.У. Математическое моделирование тлеющего газового разряда // Изв. СО АН СССР. Сер. технических наук. 1988. № 21. Вып. 6. С. 41-47.

4. Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983.

5. Григорьев Ю.Н., Горобчук А.Г. Особенности интенсификации травления кремния в CF4/O2 // Микроэлектроника. 2007. Т. 36. № 4. С. 283-294.

6. Физические величины: справочник / под ред. И.С. Григорьева, Е.З. Мейлихова. М.: Энергоатомиздат, 1991.

7. Cherrington B.E. Gaseous electronics and gas lasers. Oxford: Pergamon Press, 1980.

Статья поступила 17.06.2013 г.

Gorobchuk A. G. ON A NUMERICAL SCHEME OF EXPONENTIAL FITTING FOR SOLVING RADIO-FREQUENCY DISCHARGE EQUATIONS IN THE HYDRODYNAMIC APPROXIMATION. The numerical modeling of the RF discharge in the hydrodynamic approximation is widely used in studying processes of plasma-chemical etching of semiconductor materials. The integration of RF discharge equations with a small parameter at the high-order derivative is a difficult problem connected with the formation of boundary layers with big gradients of required values. In addition, the equations are described by a high nonlinearity and strong interrelation, which needs using special highly stable numerical methods. The basic methodological approach in developing absolutely monotonous finite-difference schemes is their regularization a special case of which is the method of exponential fitting used in this work. The proposed finite-difference scheme is considered by an example of continuity equations. The integration by the

control volume of differential expressions for particle flows yields their fmite-difference approximations of convective and diffusion components simultaneously, which provides their stable calculation at big gradients of the potential. The expressions for the flows received for the continuity equations are generalized for the electron energy equation. The introduced implicit exponential fmite-difference scheme guarantees to get the solution for big Peclet numbers with keeping positive values of electron temperature and concentrations of plasma components.Key words: radio-frequency discharge, method of exponential fitting.

Keywords: radio-frequency discharge, method of exponential fitting

GOROBCHUKAleksey Gennad’evich (Candidate of Physics and Mathematics, Novosibirsk State University, Institute of Computational Technologies of the Siberian Branch of the Russian Academy of Sciences, Novosibirsk, Russian Federation)

E-mail: alg@eml.ru

REFERENCES

1. Graves D.B., Jensen K.F. A continuum model of DC and RF discharges (1986) IEEE Transactions on plasma science. V. PS-14, pp. 78-91.

2. Lymberopoulos D.P., Economou D.J. Fluid simulation of glow discharges: effect of metastable atoms in argon (1993) Journal of Applied Physics. V. 73, pp. 3668-3679.

3. Gadiyak G.V., Shveygert V.A., Uuemaa O.U. Matematicheskoe modelirovanie tleyushchego gazovogo razryada (1988) Izvestiya SO AN SSSR: seriya tekhnicheskikh nauk. No. 21. Vyp. 6, pp. 41 -47. (in Russian)

4. Dulan E., Miller Dzh., Shilders U. Ravnomernye chislennye metody resheniya zadach s pogranichnym sloem. Moscow, Mir Publ., 1983. (in Russian)

5. Grigor'ev Yu.N., Gorobchuk A.G. Osobennosti intensifikatsii travleniya kremniya v CF4/O2 (2007) Mikroelektronika. V. 36. No. 4, pp. 283-294. (in Russian)

6. Fizicheskie velichiny. Spravochnik. Pod red. I.S. Grigor'eva, E.Z. Meylikhova. Moscov, Ener-goatomizdat Publ., 1991. (in Russian)

7. Cherrington B.E. Gaseous Electronics and Gas Lasers. Oxford, Pergamon Press, 1980.