On a numerical scheme of exponential fitting for solving radiofrequency discharge equations in the hydrodynamic approximation.pdf Исследование физики ВЧ-разряда является сложной многопараметрической задачей, которая носит преимущественно исследовательский характер. В диапазоне рабочего давления 0,1-3,0 торр относительно полное представление о распределении концентраций электронов и ионов, а также температуры электронов дает моделирование ВЧ-разряда в гидродинамическом приближении [1, 2]. Умеренная вычислительная трудоемкость гидродинамического подхода позволяет применить его для адекватного исследования процессов, реализуемых при плаз-мохимическом травлении полупроводниковых материалов. На стадии горения ВЧ-разряда, характеризующейся образованием тонких при-электродных слоев и квазинейтрального столба, из-за сильных градиентов параметров плазмы и большой проводимости проведение расчетов с применением явных численных схем для решения уравнений непрерывности электронов и ионов становится крайне неэкономичным [3]. Из-за ограничений, накладываемых условием Куранта для электронов Ke = veт / h < 1 (где ve - дрейфовая скорость электронов, т, h - шаги по времени и по пространству), шаг по времени оказывается много меньше, чем характерное время развития разряда, что может потребовать более 104 - 105 шагов по времени. Здесь эффективно применение неявных конечно-разностных схем, позволяющих использовать шаг по времени в 10 - 102 раз больший, чем дает условие Куранта. Однако в численных схемах, в которых распределение потенциала находится непосредственно из уравнения Пуассона, возникает дополнительное ограничение вида т< 1 / 4пст, где ст - проводимость плазмы [3]. Данное ограничение дает значение на порядок больше, чем условие Куранта. Кроме того, основные уравнения обычно решаются c использованием дополнительных итераций самосогласования. Само решение системы нелинейных разностных уравнений также представляет большую трудность. В настоящее время при моделировании газовой плазмы напряженность электрического поля находится из уравнения сохранения электрического тока. Для решения данного уравнения применяются консервативные по заряду разностные схемы, в которых шаг по времени ограничивается периодом лэнгмюровских колебаний ионов. При использовании чисто неявных схем в уравнении тока решение разностных уравнений требует дополнительных итераций по самосогласованию, что снижает эффективность вычислительного алгоритма. Применение в нем явных численных схем приводит к известному ограничению на шаг по времени. Поэтому для построения эффективных численных алгоритмов используют преимущества каждой из схем, которые, в конечном итоге, также обладают определенными недостатками. Интегрирование уравнений высокочастотного разряда в гидродинамическом приближении с малым параметром при старшей производной является трудной вычислительной задачей, связанной с образованием пограничных слоев с большими градиентами искомых величин, ведущих к неравномерной сходимости решения. В работе для построения безусловно монотонной разностной схемы использовался метод экспоненциальной подгонки для численного решения уравнений непрерывности электронов и ионов, а также уравнения электронной энергии, обобщенный на двумерный случай. Физико-математическая модель Рассматривается аксиально-симметричный высокочастотный разряд в гидродинамическом приближении в реакторе плазмохимического травления. Начально-краевая задача ставится в цилиндрической области, где электродами являются торцевые поверхности реакционной камеры. В расчетах используется гидродинамическая модель аксиально-симметричного ВЧ-разряда, включающая решение уравнений непрерывности для электронов и положительных ионов, уравнение баланса энергии электронов и уравнение Пуассона для электрического потенциала. Введем переменные: n ~ Vi a Di , Pi _-, i = -, Di _ -,1 _e p; ne0 Vi 0 Di 0 Ф n Te „ „ r z L ф0 Te0 rt L rt Здесь индексы e, p обозначают электроны и ионы соответственно; щ - плотность частиц; ц, Dt - подвижность и коэффициент диффузии; Te - электронная температура; ф - потенциал; f - частота активации; t - время. Плотность, подвижность, коэффициент диффузии, потенциал и температура обезразмеривались на среднюю плотность электронов ne0, характерные значения ц0, Dt0, потенциал ВЧ-электрода Ф0 и температуру Te0, соответствующую тепловой энергии электронов. Координаты нормированы на радиус цилиндра rt и межэлектродное расстояние L соответственно. Концентрации плазменных компонент находятся из уравнений непрерывности для электронов и ионов: dpi , . j J'OA_«c Т_ fL2 п _Г 1 i _ e Pi Т-дТ- + A дт - + - +--- _ PiSi, , pi ' , (1) dZ i 1 De0 Fi Щ / Dp> i _ P. Потоки частиц определяются формулами _ дф dD i pi дф dD i pi j, _-siPei|vip^----, j, _-siPei|vip^---- J% i ir7Ki д% iz i dZ dZ l = e, l = P. 0 D, п Pe, = Нижние индексы | и С относятся к направлениям радиальной и аксиальной координат. Число Пекле характеризует соотношение процессов переноса и диффузии. Источниковый член включает ионизацию исходных молекул и процесс прилипания электронов к ионам: {Pe (Da r - ra X { PeDa ri, l = e, l = P, r = exp kTe, L2 nk,r -exp 1 -J- Q. A kTef De, L2enel еФп S, = Da =- De, e0 e0 V n"Ie0, где ri - скорость ионизации, ra - скорость прилипания, Ei - энергия ионизации, k - постоянная Больцмана, va - частота прилипания, n - объемная плотность газа, ki0 - коэффициент скорости ионизации. Число Дамкелера пропорционально отношению интенсивности ионизации к диффузии. Распределение электрического потенциала определяется из уравнения Пуассона: ,2 1 д (^дф A +|ф = Г(Р.-Pp), г = I V / dZ где e - элементарный заряд, е - диэлектрическая постоянная в вакууме. Введенный параметр Г характеризует напряженность электрического поля. Электронная температура находится из уравнения баланса электронной энергии: л,. , ; дф f Л дj д Q^ J 3 TdT(pe Qe ) + A 5 дт (3) - -Х\ A2je- + je- | + ©Se = 0, дс \ ^ д| ez дz| e где компоненты потока электронной энергии имеют вид = je^Qe, jQ = je^Qe. Параметр x характеризует отношение энергии электрического поля и тепловой энергии электронов. Другой параметр © отвечает за потери энергии в реакциях ударной ионизации. Их выражения имеют вид h eФn X =' he, К, he0 = 2 kTe0 , e0 e0 где he0 - тепловая энергия электронов, h - потери энергии при ионизации. На границах расчетной области для уравнений (1) - (3) ставятся следующие краевые условия: - на нижнем электроде (0 С = 0) : Pe = 0, дD p P p = 0, ф = 0, Qe =Qe. дС - на верхнем электроде (0 , С = 1): дD p P p = 0, ф = sin(2m), Qe = Qe pZ ' дС Pe = -Yj - на оси симметрии (£ _ 0, 0 < Z < 1): _0, ^_0, дф_0, _0; д% - на боковой поверхности (£ _ 1, 0 < Z < 1): pe _ 0, ^ _ 0, 0, Q ; > e ' ^^ " --vt- " e es ~ - на перфорированной границе откачки < \ < 1, Z _ 0): д£) p дф pe _ 0, _ 0, 0, Q ; - на верхнем электроде (§o < | < 1, Z _ 1): д£) p p p дф pe _0, -^_0, ^_0, Qe _Qes; e 'z 'z es где _ r0 / rt - радиус электродов, у - коэффициент вторичной электронной эмиссии, Qes - электронная температура на электродах. ВЧ-разряд инициируется в цилиндрическом объеме, в котором низкотемпературная плазма ограничивается непроницаемой боковой стенкой в радиальном направлении. Из-за быстрой рекомбинации электронов на проводящей поверхности верхнего электрода электронная плотность полагалась равной нулю. На нижнем электроде учитывался слабый поток вторичной электронной эмиссии. Краевые условия для ионов записывались в виде диффузионных потоков к границам плазменной области. На оси симметрии задавались условия симметрии. Потенциалы и электронная температура на электродах предполагались заданными. На боковой стенке реактора, в верхней части реакционной камеры и на перфорированной границе откачки плотность электронов полагалась равной нулю, электронная температура бралась равной температуре на электродах, а для потенциала использовались мягкие краевые условия. Конечно-разностная схема Решение уравнений (1) - (3) является трудной вычислительной задачей из-за высокой нелинейности и сильной взаимосвязи уравнений. Уравнения непрерывности (1) содержат большие числа Пекле (до 104), что требует применения специальных численных методов. При интегрировании уравнений с малым параметром при старшей производной возникают сложные вычислительные задачи, связанные с образованием пограничных слоев с большими градиентами искомых величин, ведущих к неравномерной сходимости решения. В качестве одного из решений возникающей проблемы является повышение порядка точности разностной схемы и использование подробных сеток в области сильного изменения решения. Однако эффективность подобных решений существенно снижается с увеличением жесткости уравнения. Кроме того, раздельная аппроксимация конвективных и диффузионных членов (центральными или направленными разностями против потока для приближения первой производной и центральной разностью для второй производной) практически не пригодна при построении монотонных схем. Основным методологическим подходом при построении безусловно монотонных разностных схем является их регуляризация. частным случаем которой является метод экспоненциальной подгонки или интегральных тождеств. использованный в данной работе [4]. Рассмотрим суть этого метода на примере уравнений непрерывности (1). записанных для полных потоков частиц. включающих одновременно конвективную и диффузионную части. В рассматриваемой стандартной области введем прямоугольную сетку с узлами (,.Zj). i = 1.....N; j = 1.....M. Вся область разбита на некоторое количество непересекающихся контрольных объемов. Проинтегрируем дифференциальное уравнение (1) по выделенному контрольному объему. Первоначально преобразуем выражение для компоненты потока jt к следующему виду: _L dD, р, u, д, Ji, - + D, р, = n h дф ■с ■с где ul = slPe, ---. и проинтегрируем его в пределах от , до с,м: D, Щ ^ jif 1 dD, р, I =-(u, ,f+1. -+D i Pi j ъ ^ +D i Pi = - | u,d,. ln Предполагается. что потоки частиц. скорости дрейфа и коэффициенты переноса постоянны между узлами разностной сетки: D, р, I = D, р, . D,р, I = D, р, . Щ' ' ' Щ(+1 I+1 I+1 . Таким образом. приходим к равенству = 4i+1/2 . u' = \ Ji, ъ К ■+D i р, 'i + 1 r li +1 ( \ ( Л jlK Ji,. - + D i р,' ln - ln + 1/2 W +1/2 ■ V 'i +1/ 2 / V 'i +1/ 2 где Д, +1 . После преобразования получаем выражение для потока частиц Pi+1 Di 1 exP(ui 1/2 Д,'+1/2) -Pi Di _ i + 1 i +1 i +1/2 ^'+1/2 ' ' "l. К+1/2 ,) - 1 Ji, exp(u . + 1/2 ,i+1/ Данное равенство справедливо для любого j = 1.....M . поэтому окончательное выражение имеет вид Д, ) - р, D, Ji, = -u +1/2. j Pi 1 Dl eXP(ui „, Д, +1/2) -Pl D i + 1,j i +1. j i +1/2. j 1/2 i, j i. j +1/2,j exp(u, Д, ) -1 'i+1/2, j ' +1/2 Аналогичным образом вычисляются остальные потоки j, . j, . j, . -1/2, j Z'. j + 1/2 Z', j-1/2 Интегрирование по выделенному контрольному объему дифференциальных выражений для потоков частиц дает их конечно-разностные аппроксимации конвективной и диффузионной составляющих одновременно. Первые производные в уравнениях непрерывности для электронов и ионов аппроксимируются центральными разностями: k+1 k k+1 -k+1 -k+1 -k+1 -k+1 -k+1 Pi -pi Jif Jif Jif + Jif i Jir P j ii.j rii.j + A1 +1/2.J fi-1/2.j + fi+1/2.j fi-1/2.j + ^i.j+1/2 fr. J-1/2 _p pk 4 Af, 2fi Д, где потоки частиц в радиальном и аксиальном направлениях вычислялись по формулам Л+1 f\k i k+1 \ Л+1 f\ k Pi Di exp (z -Pi Di _ i +1.j i +1.J V i +1/2.j I i.j i. fi+1/2 exp(^ )-1 k+1 k+1 t\ k „ / k+1 l k+1 t\ k zi Pi Dl exP (zi )-Pi Dl ■k+1 _ _ i.J+1/2 i.J+1 i.j +1 \ i.j + 1/2 / i.j i.J zk+1 „k+1 k „„„/ „k+1 i „k+1 nk ■k+1 _ li +1/2.j ' *i+1J *i+1.J jif, " fi +1/2.J Д. +1/2 ' ' J Jir ' Ci.j+1/2 Д Z j + 1/2 ZJ + 1/2 eXP ( C/2 )- 1 Здесь индекс k относится к моменту времени Tk, при этом временной шаг составляет величину Дт _ Tk+1 - Tk . Символы zik+1 и zik+1 даются формулами k i +1/ 2. j i. j+1/2 zk _- ePe ^k +1/2.J (CDk -(Dk ) zk _- ePe ^k.J+ 1/2 (CDk -(Dk ) Zii +1/2.J _ ^ Dk +1,J l). +1/2 _ ^ Dk +1 Ci,j). ii +1/ 2. j ii. j + 1/2 В направлении f узел i+1/2 располагается посередине между узлами i и i+1 разносной сетки, которые разделены отрезком Д^ _ f{+1 . Аналогичное расположение узлов используется в направлении Z , где узел j+1/2 располагается посередине между узлами j и j+1, разделенных отрезком Д,, 1/2 _Z j+1 -Zj . Другие обозначения имеют вид Д^ _ (fг-+1 -1)/2, Д,, _(,j+1 -Zj-i)/2. Выражения для потоков jk+1 , jk+1 в узлах сетки (i-1/2, j) и (i, j-1/2) имеют аналогич- Si-1/2. j %i. J-1/2 ное представление. Значения потенциала ф*+1 j , Ck j, cfy+1, коэффициентов пеk k k k k реноса цi , ц , Di , Di и источниковый член Si в численной r г ii +1/2.J r ii.j+ 1/2 ii+1/2.J ii.j+1/2 i схеме брались с предыдущего временного слоя Tk. Для нахождения электрического потенциала решается уравнение Пуассона с помощью неявной конечно-разностной схемы типа «крест»: ( rj ,Jk k k \ фг+1,1 -фг,1 фг,1 -фг,j-1 (к к к к \ +1 -фг,j ъ,] -1 Д, ДZ ^J+1/2 ^J-1/2 у = Г(-Pki). 2 1 +- ^ ] fi Аналогичным образом интегрируются выражения для компонент потока электронной энергии. Значения электронной температуры в узлах сетки находятся из уравнения 3 Qk+V+1 _Зк „к jf1 - jf1 jf1 + jf1 3 ^eiJPei.j + a2 fi +1/2.J si-1/2.J + a2 ^i +1/2.J fi-1/2.J , 5 4 4 2 12 +x A2 ( j+1 Ek+1 + j+1 Ek+1 ) + 2 V +1/2. j ' +1/2. j -1/2.J -1/2. j/ +x( jk+1 Ek+1 + jk+1 Ek+1 )+©sk = o. 2 V eZi. j +1/2 J + 1/2 -/eZi. j-1/2 J-1/2/ e'. j где компоненты потока электронной энергии в радиальном и аксиальном направлениях вычисляются по формулам k zk+i pk+1 D exp)-sk+lpk:lDke ■k+1 _ i+1/2.j i +1.j i +1.j i +1.j \ i+1/2,j / i.j i.j i.j Jk _ л " Zi +1/2, j Д " +1/2 Zi+1/2 exp (zk+1 )-1 V i +1/2. j / zk+1 sk+1 pk+1 % exp (zk+1 )-»k+1pk+1 D i.j+1 i.j+1 \ i.j+1/2 / i.j i.j i jk+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 J\i,j+1/2 " Д^ /2 4+1/2 exp (z^ 1/2)-1 Компоненты электрического поля определяются как t?k+1 /_k+1 k+1 \ it t?k+1 ik+1 k+1 \ < л +1/2,j ="(^+W ) 7 Д^,+1/2. EC,j+1/2 =-(^k J+1 ) / ^^ Выражения для компонент потока электронной энергии jk1 . jk1 и элек- Z'-1/2, j Ci, j-1/2 трического поля Ek+1 . E^+1 имеют аналогичное представление. Вследствие i-1/ 2. j i. j-1/ 2 локального выполнения законов сохранения для конечных контрольных объемов решение даже на грубой сетке удовлетворяет точным интегральным соотношениям. Предложенная численная схема сохраняет положительные значения концентраций плазменных компонент и электронной энергии. Учитывая наличие больших градиентов электрического поля. для устойчивого расчета потоков электронов и ионов на границах области (на электродах и боковой стенке камеры) использовалась формула линейной интерполяции. позволяющая выразить поток на границе через значения потока внутри области. рассчитанного по экспоненциальной схеме. например. на верхнем электроде: jk+1 _ jk+1 Дг Jkr Jkr jk+1 = jk+1 + Z j +1/2 Zi, j +1/2 Zi, j-1/2 JkZi.j+1 " JkZi.j+1/2 2 ' Дг • z j Каждое из уравнений в сеточной области записывается в виде системы линейных алгебраических уравнений. имеющей матрицу большой размерность с разреженно-упорядоченной структурой. Учитывая небольшое количество переменных в строке или столбце. для решения таких СЛАУ эффективно использовать классические итерационные или проекционные методы (с предобуславливанием) на подпространствах Крылова. В работе СЛАУ решались итерационным методом Гаусса - Зейделя, показавшим хорошую эффективность. Результаты и обсуждение В расчетах моделировался ВЧ-разряд в плазмохимическом реакторе травления радиальной схемы при рабочем давлении p = 0.5 торр [5]. Радиус реакционной камеры выбирался равным rt = 30 см. а межэлектродное расстояние L = 3.5 см. ■ к+1 k+1 J»_ , Jar + ^i. j+1/2 Si.j-: 4 ВЧ-разряд инициировался на рабочей частоте f = 13,56 МГц при напряжении на электродах Ф0 = 110 В. Другие параметры выбирались следующими: плотность газа n = 3-1016 см-3, средняя плотность электронов ne0 = 6-109 см-3, электронная температура на электродах и боковой стенке Tes = 0,5 эВ. Другие размерные величины, входящие в уравнения, брались из работ [1, 2, 5, 6]. Начальные распределения концентраций электронов и ионов выбирались в виде аналитического решения, полученного в диффузионно-дрейфовом приближении для ВЧ-разряда в цилиндрической области [7]: pe0 =р 0 = 3,63796 J0 (2,405^) sin(nZ). Начальное распределение электронной температуры полагалось постоянным (рис. 1). Рис. 1. Изолинии электронной плотности ne -10 10 см 3 в плазмохимическом реакторе Моделирование ВЧ-разряда с гармоническим током разряда требует решения нестационарных уравнений переноса зарядов и баланса энергии электронов. В расчетах рассматривалось до 104 периодов ВЧ-разряда, когда колебания электронной плотности и температуры выходили на периодический режим. Из-за высокой нелинейности исходных уравнений, включающих решение уравнения Пуассона, полученные неявные численные схемы являются условно устойчивыми. В этой связи шаг по времени выбирался в диапазоне Дт « 5 -10-6 -10-4 , при этом размер сетки составлял 40x20 узлов по пространственным переменным. Тестирование алгоритма проводилось на основе сравнения численного решения с аналитическим, полученным в диффузионно-дрейфовом приближении для ВЧ-разряда между двух поверхностей [7]: pe0 = 1,57079sin(nZ). При этом на верхнем электроде (0

Скачать электронную версию публикации