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

Размерные модели реакторов ГХОГН. 3-D(x,y,z) модельРазработанная 3-D(x,y,z) модель включала в себя взаимосвязанные транспортный блок,блок газофазной химической кинетики, блок процессов на поверхности нити и подложки. Длятечений вязкого, сжимаемого и реагирующего газа в прямоугольной системе координат (x,y,z)численно решалась система уравнений сохранения массы газа, компонент смеси, импульсов иэнергии:3-D(x,y,z) модель уравненияr∇t ρ = −∇( ρv )r∇t (ρu) = −∇(ρuv) − ∇z p + ∇Wz + ρg(1)r∇t (ρv) = −∇(ρvv) − ∇x p + ∇Wx(3)r∇t (ρw) = −∇(ρwv) − ∇y p + ∇Wy(4) kr∇t (ρε ) = −∇(ρεv) − p(∇ zu + ∇ x v + ∇ y w) − ∇Qλ − ∇ ∑ hi jiD  + Wε i=1r∇ t (ni ) = −∇ ni v + jiD / mi + Si − Li ni()(2)(5)(6)42Для замыкания этой системы использовались уравнения состояния идеального газа,термическоеp=N×k×T(7)и калорическое (для газовой смеси из nk компонент 1≤i≤nk с их энтальпиями образования ∆Hf0iи температурно-зависимыми теплоемкостями CVi(T) [163])ρε = ∑ ρiε i = ∑ kni (∆H f 0 + TCVi ) / Rii(8)iгде R – универсальная газовая постоянная (R=1.987262 кал/(моль К) для ∆Hf0i в кал/моль),k=1.38×10-16 эрг/К - постоянная Больцмана (R[эрг/(моль К)]=kNA, NA–число Авогадро).

Вуравнениях (1-6) частная производная∇b =∂,∂qq = t, z, x, y; дивергенцияr∇(ρv ) = ∂ (ρu ) ∂z + ∂ (ρv ) ∂x + ∂ (ρw) ∂y , ρ - плотность газа, u, v, w - компоненты скорости внаправлении z, x, y , p - давление, Wz, Wx, Wy и Wε – вязкостные члены в уравнениях импульсови энергии[159,160],24Wz =  μ∇ z u - μ ⋅ (∇ x v + ∇ y w), μ ⋅ (∇ xu + ∇ z v ), μ ⋅ (∇ y u + ∇ z w)3342Wx =  μ ⋅ (∇ x u + ∇ z v ), μ∇ x v - μ ⋅ (∇ z u + ∇ y w), μ ⋅ (∇ y v + ∇ x w)3342Wy =  μ ⋅ (∇ y u + ∇ z w), μ ⋅ (∇ y v + ∇ x w), μ∇ y w - μ ⋅ (∇ x v + ∇ z u) 33µ = µ(T) - коэффициент вязкости смеси. Для возможности увеличения шага по временичисленного метода на начальных этапах счета вводилась исскуственная вязкость (увеличеннаяна порядок µ).

g – ускорение силы тяжести (сила тяжести –рассматриваемых реакторах [17]),ε[эрг/г]малозначительный член в- удельная энергия газовой смеси, εi и hi -удельные энергия и энтальпия i-ой компоненты. Потоки тепла Qλ = (− λ∇ zT ,−λ∇ xT ,−λ∇ yT )T(9)по направлениям вдоль осей z, x y, λ=λ(T) - коэффициент теплопроводности, T - температурагаза, ni и mi - концентрация и масса i-ой компоненты. Диффузионный поток i-ой компоненты восновном газе k [161]ji = − mi Di ⋅ N ⋅ (∇X i +DkTi⋅ ∇T )T(10)Xi=ni/N, N - концентрация газа, Di= (mk/(miXi + mkXk))Dik - эффективный коэффициент диффузии,Dik [см2 с-1] и kTi- коэффициент бинарной диффузии и термодиффузионное отношение[161,162]. Бинарные коэффициенты диффузии i-ой компоненты в доминирующей (k-ой)43компоненте (Н2 в рассматриваемых реакторных смесях) вычислялись по формулам Dik =ai∙Tb/p,b~1.5-1.7, детали вычислений Dik приведены в главе 6.

Si и Lini - скорости образования и гибелиi-ой компоненты в химических реакциях.Система газофазных химических реакций обрабатывалась с помощью разработанноготранслятора химических реакций, который автоматически формировал правые части (Si - Lini) иЯкобиан системы уравнений химической кинетики. Для H/C смеси, для которой использовалсямеханизм GRI Mech 3.0 [163], в базовый вариант входило 39 обратимых реакций для 14компонент: H, H2, CHx, x=0-4 (для метилена, кроме основного триплетного состояния, отдельнорассматривался низколежащий синглетный уровень СH2(S)) и C2Hy, y=1-6.

Кроме того,использовались различные расширения химической кинетики с добавлением реакций для С2(a)и С2(X) компонент, азотных компонент для Н/C/N смесей (§2.4), бор-содержащих компонентдля H/C/B/O смесей (§2.5). Константы скоростей обратных реакций определялись также спомощью транслятора химических реакций по термодинамическим и термохимическимданным [163,164].Системы уравнений химической кинетики зачастую являются так называемымижесткими системами, в которых малые по величине, но быстро меняющиеся компоненты могутсущественно ограничивать допустимую величину временного шага явной схемы. Для решенияуравнений химической кинетики применялись и тестировались различные специальныеметоды: неявный метод Гира [165], разработанный метод автоматического выделения«быстрых» компонент [166].

Суть последнего метода заключается в следующем. На каждомвременном шаге используется неявная схема Эйлера. Получающаяся при этом системанелинейных уравнений решается итерационным методом. Сначала делается несколькоитераций Зейделя, затем только для тех компонент, для которых не была достигнутасходимость, делается одна - две итерации Ньютона с аналитическим вычислением Якобианаподсистемы меньшей размерности. Эти методы позволяли успешно решать отдельно системуобыкновенных дифференциальных уравнений химической кинетики (0-D модель без всякихвнешних возмущений – например, вклада fT транспортных процессов), но оказались не оченьпригодны для использования в 2-D и 3-D моделях из-за проблем с консервативностью ибольшими вычислительными затратами.

Поэтому был разработан простой практический метод,пригодный для встраивания в полноразмерные модели. А именно, явная схема для большинствакомпонент(«медленных»компонент)иполунеявнаядлябыстрых,какправило,малонаселенных, компонент:ni(t+dt)= ni + dt∙(Si - Lini + fT) ,если dt∙Li<Acni(t+dt)= (ni + dt∙(Si + fT))/(1 + dt∙Li) , если dt∙Li≥Ac(«медленные» компоненты)(«быстрые» компоненты)Величина порога Ac, разделяющего «медленные» и «быстрые» компоненты, варьировалась и44после серии тестовых расчетов была выбрана фиксированнойAc=0.05 для всехрассматриваемых задач.Нестационарные уравнения (1)-(6) с учетом термического и калорического уравнений(7,8) состояния идеального газа, граничных и начальных условий решались численно доустановления стационарного режима.

Вместе с полунеявной схемой для компонент смеси врасчетахуравненийсохраненияимпульсовиэнергиииспользоваласьявнаясхемаконсервативного конечно-разностного метода, подобная схеме метода крупных частиц [167] сострогим соблюдением балансов импульсов и энергии в расчетной ячейке. Плотность газа,концентрации компонент смеси, внутренняя энергия, а также температура, давлениеприписывались к центрам ячеек (круным частицам), а скорости и все потоковые члены - ксоответствующим границам (центрам сторон ячейки). Размеры пространственных ячеек были,как правило, 0.1×0.1 см во всех 2-D моделях и 0.1×0.3×0.1 см в 3-D(x,y,z) моделях реакторовГХОГН.

Увеличивать размеры ячеек в перпендикулярных оси нити направлениях (x,z,y – вдекартовых координатах и r – в цилиндрических 2-D(r,z) координатах) не позволяютзначительные градиенты температуры газа Т и концентраций компонент в этом направлении.Ось z декартовой системы координат во всех моделях была перпендикулярна подложке и нити(оси спирали) и направлена от центра масс спирали (системы нитей, и тогда здесь z=0) к центруподложки (или наоборот, от центра подложки, которому тогда и будет соответствовать z=0), осьх направлена вдоль подложки в плоскости, перпендикулярной нити (оси спирали), ось yнаправлена вдоль нити (оси спирали), как показано на схеме ГХОГН реакторов с одиночнойспиралью (рис.2.1a) и многонитевой системой параллельных нитей (рис.2.1b).Перейдем к постановке граничных условий (ГУ) и особенностям описания геометрии иприсутствия горячей нити в модельном реакторе в такой системе координат.

Для скоростей газана твердых поверхностях использовались стандартные условия непротекания и прилипания.Для модельной геометрии ввода газа в камеру (например, для 2-D(r,z) моделей центральнойтрубки на оси камеры реактора или коаксиального кольца) скорость влета исходного (прикомнатной температуре) рабочего газа в реактор бралась такой, чтобы обеспечить заданныйэкспериментальный расход газа. При постановке ГУ на подложке для компонент смесисчиталось, что в поверхностных реакциях участвуют атомарный и молекулярный водород игазофазные предшественники алмаза (ПА, компоненты, из которых строится АП). Для них вследующем подразделе 2.1.3 рассматривается кинетика газофазно-поверхностных реакций имеханизмы роста АП.Детальные параметры и конфигурация модельного реактора выбирались близкими кпараметрам экспериментальных реакторов с различными геометриями ГН: ГН в виде спирали(рис.2.1a), одиночной длинной нити или системы параллельных нитей (рис.2.1b).

Корректное45воспроизведение в модели реальных объемов реактора (удаленности холодных стенок),объемов горячей зоны (например, цилиндрической области, внутри которой заключенаспиралевидная ГН) и теплопроводностных потоков с поверхности ГН имеет принципиальноезначение для расчета реальных профилей температур (и, значит, концентраций компонентсмеси) в реакторах ГХОГН. В этом смысле для каждого моделируемого реактора важен выборадекватной размерной модели 2-D(x,z) (двумерной декартовой для длинных нитей иподложкодержателей, когда их длины вдоль оси y много больше, чем характерные размеры подругим осям, ширина подложки по оси х, зазор Lfs, и т.д.), 2-D(r,z) (цилиндрической модели длякольцевой ГН или спирали ГН в отсутствии подложки, уравнения в цилиндрической геометрииприведены в главах 3 - 5) и трехмерной 3-D(x,y,z) модели, позволяющей наиболее точновоспроизводить реальную геометрию любого реактора, но наиболее требовательной в смыслевремени счета и ресурсов компьютера.