Трехмерное численное моделирование процесса переноса паров кремния в условиях среднего вакуума в задаче с цилиндрической симметрией

В последние десятилетия использование композиционных материалов становится все более технологичным и разнообразным. В результате их совершенствование и производство бурно развиваются во всем мире [1 –5] . Можно предполагать, что композиционные материалы будут находить все новые области применения. Одним из активно используемых в технике композитов является насыщенное кремнием углеродное волокно, обладающее высокими антиокислительными, прочностными и герметизирующими свойствами [6 –11] . Такой композит получается путем жидкофазного силицирования [12, 13] . Однако при реализации данной технологии присутствующие в системе капиллярные силы не позволяют проводить глубокую и равномерную пропитку всего объема материала, необходимую для достижения требуемых параметров проницаемости.

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

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

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

Статья опубликована в открытом доступе по лицензии CC BY 4.0

технологического процесса время порядка нескольких часов и не требует приближения тиглей к изделию на максимально близкое расстояние. Ранее авторами настоящей работы предложено как объяснение перечисленным выводам [19] , так и более сложная физико-математическая модель процесса переноса паров кремния в разреженной среде инертного газа (как правило, аргона), в которую добавлен еще один механизм переноса паров кремния. Этот дополнительный механизм транспорта имеет конвективную природу, и его учет позволил прийти к новому уравнению нелинейной диффузии, которое сначала было отработано на системах с более простой геометрией. С использованием этого уравнения были получены точные стационарные решения одномерных распределений паров кремния в задачах с декартовой и цилиндрической симметрией [20] , решены двумерные задачи с модельными вариантами расположения зеркала расплава и изделия [21] . Была решена задача в трехмерной постановке для тигля конечных размеров и реторты, имеющей форму прямоугольного параллелепипеда [22] . В работе [23] проанализирован случай неоднородного распределения температуры в реторте. Показано, что это привносит дополнительный вклад в транспортные свойства паров кремния в зависимости от величины и направления градиента температуры в системе. Полученные результаты позволили сформировать общую непротиворечивую картину динамики концентрационного фронта паров кремния и стационарных распределений полей их концентрации в реторте в ходе силицирования для простейших конфигураций системы.

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

• 2.    Постановка задачи

Ранее в работах [19 –22] перенос паров кремния в разреженной среде аргона с учетом конвективного механизма транспорта описывается модифицированным уравнением:

∂C ∂t

- D c ( ∇ C) ² =D∆C.

Это квазилинейное дифференциальное уравнение в частных производных второго порядка параболического типа [24], в котором C — массовая доля (или массовая концентрация) газообразного кремния [25], D — коэффициент диффузии [м2/с], Dc — параметр конвективного переноса, имеющий аналогичную размерность. Слагаемое, квадратичное по градиенту концентрации, не вводится в уравнение искусственно, а естественно возникает из выражения (u∇) C, отвечающего за конвективный перенос в системе. Следовательно, параметры D и Dc математически не связаны, и каждый из них определяется собственными физическими факторами [19–22]. В эксперименте давление стороннего газа аргона в реторте соответствует условиям среднего вакуума и составляет примерно pa = 100 Па. Давление паров кремния не может превышать значения насыщения, что в рассматриваемых условиях при температуре порядка 1800–1900 К составляет примерно psi = 10 Па [26–28]. Таким образом, в отношении парциальных давлений компонентов смеси имеет место достаточно сильное условие: pa ≫psi .

Параметр конвективного переноса в уравнении (1) , появляющийся при выводе уравнения, определяется как

(см. [20 –22] ):

_{D =}3 kT

8 σ12p

πkT 2m12

3 (kT) ^{3 / 2}

8 d ² p ^√ πm0

=0.7м ² /с,

1   2 ^{3 / 2 √}m0kT

η=

3 π        d ²

=1.1 · 10 ^-4 м ² /с.

Здесь σ12 — эффективное сечение рассеяния для двух частиц; m12 — приведенная масса атомов кремния и аргона; k — постоянная Больцмана; p — давление аргона. Для двух примерно одинаковых по массе и размерам частиц имеем: σ12 = πd ² , где m12 = m0 /2, m0 — масса одного атома кремния ( ∼ 4.7 · 10 ^-26 кг), d— диаметр атома аргона (d= 2.3 · 10 ^-10 м). Величины получены для давления аргона p a = 100 Па и температуры T = 1800 К, однако вариация этих параметров в рабочих интервалах не изменяет на порядки коэффициенты диффузии и динамической вязкости.

Согласно уравнению Менделеева–Клапейрона, плотность насыщенного пара, находящегося над зеркалом расплава, является малой величиной:

ρ = ^{p si µ si} = 1.87 · 10 ^-5 кг/м ³ .

si RT

Оценка плотности кремния выполнена также при T = 1800 К и давлении насыщения p si = 10 Па.

Определим еще один важный параметр — пористость аргона ϕ. Заметим, что она вычисляется достаточно точно, без каких-то излишних и необоснованных предположений, и равна доле пустого пространства в объеме, занимаемом несущей средой. Пусть v ₀ — объем одного атома аргона, n a = p a /kT — число частиц аргона в единице объема расчетной области (парциальное давление аргона p a = 100 Па из эксперимента). Тогда имеем значение пористости, стремящееся к единице: ф = 1 — v o n a = 0.999999997. Вследствие этого в условиях среднего вакуума атомы аргона выступают в роли рассеивающих центров и занимают чрезвычайно малый объем. Проницаемость пористой среды для газовой смеси аргона и паров кремния вычислим по формуле Козени–Кармана, что обоснованно, так как она выводится из общих гидродинамических соображений и применима для сред с произвольной пористостью [23] :

κ 0 ϕ ³

К =------д,

(1 — ф) ²

где константа Козени–Кармана κ ₀ в формуле пропорциональна квадрату характерного размера обтекаемых препятствий, роль которых играют атомы аргона: к ₀ = a ² ; a = 1.4 • 10 ^-10 м — характерный размер атома. Таким образом, проницаемость оказывается равной: к =2 ^ 10 ^-3 м ² .

Подстановка всех полученных значений в формулу (2) для параметра конвективного переноса дает оценку: D c = 54 м ² /с.

• 3.    Стационарное решение в цилиндрической системе координат

  

  Рис. 1. Геометрия модельной задачи и система координат

  
  

  Для понимания сути проблемы сначала рассмотрим модельную одномерную задачу парофазного силицирования изделия в цилиндрической постановке. Пусть внешняя оболочка радиусом R ₂ является зеркалом расплава в реторте и генерирует пары кремния, как представлено на рисунке 1. Легко показать, что действие силы тяжести не оказывает существенного влияния на перенос паров кремния в реторте. Как сказано выше, плотность разреженного газа мала и вычисляется через уравнение состояния идеального газа. Для парциального давления 10 Па (давления насыщения парами кремния при температуре 1800 К) плотность составляет около p si = 1.87 • 10 ^-5 кг/м ³ . Объемные плотности сил давления в печи размерами L порядка метра и силы тяжести с вектором ускорения g ( | g | = 9.81 м/с ² ), соответственно, равны:

I^P si |~ | p si | И 10 Па/м, | p si g | И 1.83 • 10 ^-4 Па/м.

Таким образом, имеем: |∇ p _si | ≫ | ρ _si g | . В рассматриваемом диапазоне температуры влиянием ее неоднородности в реторте также можно пренебречь [23] .

Цилиндр радиуса R ₁ (Рис. 1) представляет собой изделие, изготовленное из пористого материала, полностью поглощающее пары кремния. Пусть на поверхности расплава поддерживается значение концентрации C s , соответствующее состоянию насыщения [30] . При давлении насыщенного пара кремния в рассматриваемых термодинамических условиях и соотношении 1:10 парциальных давлений кремния и аргона значение концентрации насыщения примем равным C _s = 0.08. Сначала решим задачу в предположении однородности по аксиальной координате z и отсутствия зависимости от угла α .

Итак, требуется найти решение одномерного стационарного уравнения (1) , которое в цилиндрических координатах с учетом симметрии (d/dz = 0 и д/да = 0) принимает вид:

-

dC \2 dr

ψ d r dr

при следующих граничных условиях

r = R 1 : C = 0;     r = R 2 : C = C s .

В уравнение (3) введено обозначение ф = D/Dc — безразмерный параметр задачи, показывающий отношение диффузионного и конвективного вкладов в массоперенос.

Интегрирование краевой задачи (3) и (4) дает специфическое суперлогарифмическое решение вида [21] :

C (r)= ф1п

exp(C s /ф)1n(r/R1) — 1n(r/R2) 1n(R2/Ri)

В предельном случае Dc = 0 уравнение (1) с теми же граничными условиями приводит к более простому логарифмическому решению:

C (r) = C s

ln(r/R 1 ) ln(R ₂ /Ri).

Табуляция решений (5) и (6) представлена на рисунке 2 (кривые 4 и 3 соответственно) в зависимости от относительной радиальной координаты (r — R 1 ) / (R 2 -R 1 ). При построении графиков приняты следующие значения параметров: ^ = 0.016, C s = 0.08, R 1 = 0.2 ми R 2 = 0.9 м (для решений 1 и 2 в прямоугольной системе координат считалось, что горизонтальная координата x = r). Видно, что решение (5) , как и в случае плоской геометрии [19] , представляет собой сильно выпуклую функцию, которая описывает распределение концентрации паров с резким градиентом при r = R1, то есть на входе в изделие. Оценки показывают, что столь выраженная динамика изменения концентрации обеспечивает высокую скорость массопереноса кремния, достаточную для насыщения пор углеродной матрицы за разумное время, хорошо согласующееся с экспериментальными данными. Эта зависимость значительно отличается по своей структуре от аналогичного решения, полученного ранее на основе чисто диффузионной модели [17, 18] . Упрощенная модель [17, 18] в прямоугольной системе координат приводит к линейной зависимости концентрации от координаты (линия 1 , Рис. 2) , дает на порядок меньшие значения массопотока и несостоятельна в плане предсказания времени силицирования углеродных преформ.

( r - R 1 ) / ( R 2 - R 1 )

Рис. 2. Кривая 1 – линейное решение (классическая диффузия в прямоугольной системе координат); 2 – решение при условии D c = D (эффект нелинейной диффузии присутствует, но он мал и по величине такой же, как в типичном квазилинейном уравнении, то есть отличие вкладов линейного и квадратичного слагаемых в уравнении (1) не учитывается; также не принимается во внимание топологический эффект схождения потока к оси симметрии); 3 – решение (6) при D c =0 (эффект нелинейной диффузии отсутствует, учитывается только топологический эффект схождения потока к оси симметрии); 4 – решение уравнения (1) при D c ≫ D (принимается во внимание как реалистичный вклад нелинейной диффузии, так и топологический эффект схождения потока к оси при наличии цилиндрической симметрии)

Как показывает рисунок 2, для системы с цилиндрической симметрией отличие между модифицированным и классическим уравнением диффузии (кривые 4 и 3 соответственно) становится менее очевидным, и их сравнение, особенно при трехмерной постановке, возможно лишь количественно. Однако видно, что конвективный механизм привносит существенный вклад в массоперенос кремния и обеспечивает на входе в изделие производную, практически на порядок превышающую значение для чисто диффузионного транспорта.

Для предельного случая D c = D уравнение (1) известно достаточно давно; иногда его называют уравнением Л.С. Лейбензона [31] . Впервые необходимость его использования возникла при описании процессов фильтрации в пористых средах, и оно фактически вытекает из закона сохранения массы. Применительно к задаче о распространении тепла при теплопроводности, зависящей от температуры, подобное квазилинейное уравнение рассмотрено в [32] . В общем контексте (с учетом более высокой степени нелинейности) без привязки к конкретным граничным условиям его математические свойства анализировались в [33 –35] . Однако сейчас нелинейное слагаемое в уравнении (1) не является малым. Оно имеет не вторичную по отношению к диффузии, а самостоятельную физическую природу. Как уже отмечено, данное слагаемое моделирует дополнительный конвективный массоперенос в системе, связано со слагаемым ( u V ) C и превалирует над диффузионным массопереносом. Параметры D и D c не зависят друг от друга, при этом, с точки зрения неравновесной термодинамики, градиент концентрации кремния можно интерпретировать как дополнительную термодинамическую силу, вызывающую поток его паров. В случае сложной трехмерной задачи с большим числом источников, которые описываются кусочными граничными условиями, найти точное решение уравнения (1) не представляется возможным, поэтому требуется численное моделирование.

Для математического представления процесса переноса кремния при парофазном силицировании, конечно, можно сформулировать более полную постановку задачи [36] , которая будет учитывать эффект оттеснения аргона парами кремния. Работа [36] содержит полную газодинамическую задачу перераспределения паров кремния, когда замыкание производится с помощью уравнения Навье–Стокса. Показывается, что оттеснение аргона имеет место, однако отличие концентраций аргона вблизи зеркала расплава и на пористом образце составляет всего 5%. Этот факт подтверждает применимость предлагаемого далее подхода.

• 4.    Трехмерная задача

Рассмотрим далее более реалистичную постановку задачи, в которой изделие и реторта имеют цилиндрическую форму, а на дне реторты находится один или несколько тиглей конечных размеров (Рис. 3а , б ). Конфигурация с

( a )

Рис. 3. Схематическое изображение системы с одним ( а ) и двумя тиглями ( б ) в трехмерной постановке; стрелками обозначен сходящийся поток паров кремния

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

По-прежнему на зеркале расплава Γ устанавливается концентрация насыщения C s . На поверхности изделия выполняется условие полного поглощения. На вертикальной наружной стенке реторты, верхней границе и дне ставится условие непроницаемости — отсутствие потока через эти границы. В силу симметричного расположения тиглей вокруг изделия в опыте достаточно провести моделирование в сегменте с учетом периодичности решения (в том числе с попарной установкой тиглей) по азимутальной координате.

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

r = R 2 : ^d C = 0;

∂r

r = R i : C = 0;

ас z = °ИГ 1:^ =0;

z = 0( G Г ): C = C s ;     z = H :    =0;

∂z

C (r,a,z) = C(r,a+Aa,z).

• 5.    Методика численного решения

Прямое 3D-моделирование уравнения (1) с граничными условиями (7) для поля C (r,a,z,t) и поиск стационарных решений C (r,a,z) выполнялись в авторском пакете CrystarPack [37, 38] с использованием метода конечных разностей. Пакет, написанный на языке программирования Фортран-90 и реализованный как DLL-библиотека, предназначен для решения задач механики сплошных сред и значительно сокращает время разработки расчетных программ.

Пакет CrystarPack поддерживает 2D/3D структурированные сетки с регулярным либо адаптивным шагом, разнесенные узлы сетки и метод деформированных координат для сложной геометрии. Доступны схемы конечных объемов (1–3-го порядка точности), а также конечных разностей с разностными производными произвольного порядка и точности, явные/неявные схемы и адаптивный шаг по времени. Физическая модель (уравнения и граничные условия) задается компактно через процедуры пользовательского пакета, что обеспечивает ее гибкость по сравнению с «закрытыми»коммерческими аналогами.

Вычислительным ядром является модифицированный метод Ньютона–Рафсона с механизмами глобальной сходимости и декомпозиции по полям/областям. Матрица Якоби рассчитывается численно, с экономией времени и рабочей памяти. Для решения СЛАУ доступны итерационные (GMRES+ILU по умолчанию, SPARSKIT2, LIS) и прямые решатели (MUMPS, HSL/M41). Пакет CrystarPack заполняет нишу между коммерческими программными комплексами и разработкой кода с нуля. Аппроксимация с помощью неявной схемы позволяет значительно сократить число шагов по времени (до 50–100) при реализации нестационарного расчета (например, до выхода решения на стационарный режим) по сравнению с явной схемой, требовавшей до 10 тысяч шагов. Кроме того, реализованный авторами итерационный метод Ньютона–Рафсона дает возможность напрямую решать задачу поиска стационарного режима. Таким образом, в обсуждаемой работе предлагается устойчивая и вычислительно эффективная методика расчета. Применение многопроцессорных вычислительных технологий сократит время расчета еще существеннее.

Расчеты осуществлялись на компьютере с процессором Intel® Core™ i5-4690 х 4 с 8.0 Гбайт оперативной памяти. Для аппроксимации применялась трехслойная неявная схема второго порядка точности как по времени, так и по пространству. Граничные условия дискретизировались также согласно формулам второго порядка точности. Число узлов по радиальной, азимутальной и аксиальной координатам варьировалось от 51 х 61 х 21 до 101 х 181 х 41. Вычислительные эксперименты на различных сетках показали хорошую точность и сходимость результатов.

Рис. 4. Характер концентрации насыщения на стадии установления в центре вычислительного домена в зависимости от фактора размера сетки; сегмент с одним тиглем (кривая 1 ), сегмент с двумя тиглями ( 2 )

На рисунке 4 представлены графики изменения концентрации насыщения в центре реторты при одном и двух тиглях в зависимости от фактора размера сетки N, который входит как множитель в следующие формулы вычисления числа узлов: (10N + 1) х (12N + 1) х (4N+1) для конфигурации с одним тиглем; (10N +1) х (18N +1) х (4N+1) для конфигурации с двумя тиглями.

В задаче с одним тиглем в рабочей области его размеры принимались равными 15.7 х 14 см, где первое число отражает длину тигля вдоль азимутальной координаты по внешней дуге, а второе — по радиальной координате. Зеркало расплава располагалось на дне реторты у внешней вертикальной границы (Рис. 3а ). Под тигель с расплавом кремния выделялась расчетная сетка размером 21 х 21 узел. В конфигурации с двумя находящимися рядом тиглями (Рис. 3б ) размер каждого из них был таким же(15.7 х 14 см), и аппроксимация осуществлялась с тем же числом узлов на каждый тигель.

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

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

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

Параметры расчетной области соответствовали их характерным значениям в печах для проведения процесса парофазного силицирования. Радиусы изделия и реторты составляли R ₁ = 0.2 м и R ₂ = 0.9 м, а их высота для простоты принята одинаковой: H = 40 см. Численное моделирование проводилось для одного или двух тиглей с углами раствора сегментов Да 1 = 60 ° и Да 2 = 90 ° . Очевидно, что первый вариант соответствовал шести тиглям вокруг изделия, второй — восьми. Форма тиглей имела вид трапециевидного сегмента, что обеспечивало совпадение границ с узлами сетки в цилиндрической системе координат. Уравнение нелинейной диффузии (1) решалось при значениях параметра конвективного переноса D c = 52.8 м ² /с и при D c = 0 (в предельном случае классической диффузии). Значение параметра выбиралось таким же, как в работе [22] , чтобы иметь возможность выполнить хотя бы качественное сравнение с результатами решения для плоской геометрии.

Поле концентрации в горизонтальном сечении вблизи дна (на высоте h = 1 см) представлено на рисунке 5 в момент времени t = 0.01 c. Расчетная сетка содержала 101 х 121 х 41 узлов по радиальной, азимутальной и вертикальной координатам для конфигурации с одним тиглем и 101 х 181 х 41 узлов для конфигурации с двумя тиглями. Видно, что за это небольшое время над тиглями уже сформировался концентрационный фронт паров кремния, однако он еще не успел дойти до изделия.

То, что концентрация достигает максимума над самим зеркалом расплава, ожидаемо, так как именно там происходит парообразование. Наблюдается ярко выраженная граница концентрационного фронта, что характерно для обобщенного уравнения нелинейной диффузии u t = ∆(u ^m ) при m > 1 [34] .

Любопытно, что цилиндрическая симметрия поля концентрации слабо нарушается из-за локальности источников. Диффузия даже за такой малый промежуток времени размывает неоднородности, связанные с дискретным расположением источников кремния. В результате на этом этапе в реторте уже имеет место сходящийся к оси симметрии регулярный концентрационный фронт, граница которого с хорошей степенью точности обладает цилиндрической симметрией. Ввиду сильной разреженности среды и высокой температуры (см. зависимость кинетических коэффициентов от этих параметров) концентрационный фронт быстро достигает изделия, и уже ко времени t ~ 1.0 c поле концентрации перестает меняться и приходит в состояние установления. На рисунке 6 приведены поля концентрации паров кремния при одном и двух тиглях вблизи дна. Несмотря на разное число тиглей на единицу длины дуги сегмента, профили качественно похожи. Тем не менее значения производных на входе в изделие различаются: 3.828 м ^-1 для одного тигля и 4.068 м ^-1 для двух. Во втором случае тиглям отвечает бо´ льшая общая площадь зеркала расплава, поэтому поток кремния при двух тиглях выше, чем при одном. Вполне естественно, что над тиглями наблюдается достаточно большое количество пара, и это

Рис. 5. Поле концентрации в переходный момент времени t = 0 . 01 св горизонтальном сечении на высоте h = 1 см; сегмент с одним тиглем ( а ), сегмент с двумя тиглями ( б )

( a )

Рис. 6. Поле концентрации в стадии установления в горизонтальном сечении на высоте h =1 см; сегмент с одним тиглем ( а ), сегмент с двумя тиглями ( б )

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

( a )

Рис. 7. Поле концентрации в стадии установления в горизонтальном сечении на высоте h = 20 см; сегмент с одним тиглем ( а ), сегмент с двумя тиглями( б )

Движение паров кремния фактически представляет собой локализованные над его расплавом восходящие струи. Сначала движение вверх характеризуется сильно неоднородным распределением концентрации кремния, но с течением времени профиль концентрации сглаживается и приходит к стационарному виду. Поле концентрации над тиглями в вертикальном сечении при r = 0.83 м в различные моменты времени представлено на рисунке 8.

Подобное распределение паров кремния в реторте говорит о том, что кремний в реторте есть везде, а не только над тиглями (см. Рис. 9а и Рис. 9б ). Этот кремний активно проникает в поры по всей поверхности изделия. В свою очередь это говорит о том, что не имеет смысла приближать тигли к изделию на максимально близкое расстояние. Если силицирование на каких-то участках изделия не происходит, то это уже не связано с отсутствием кремния вблизи его поверхности, а имеет совсем другое объяснение. В этой ситуации ключевую роль начинает играть

( a )

( в )

( б )

( г )

0.2 0.4 0.6 0.8 ос

Рис. 8. Поля концентрации кремния над тиглями в вертикальном сечении изделия при r = 0 . 83 м в разные моменты времени t , c: 0.1 ( а ), 0.2 ( б ), 0.5 ( в ), 1.0 ( г )

( a )

Рис. 9. Вертикальное распределение концентрации кремния над тиглем ( а ) и в пространстве между тиглями ( б ) в разные моменты на отрезке времени 0 ⩽ t ⩽ 1 . 0 c при шаге ∆t = 0 . 2 c (кривые 1 – 6 )

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

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

Таблица. Значения производных на поверхности изделия вдоль вертикали в единицах 1/м

Высота, см Число тиглей (вид диффузии)	1	10	15	20	25	30
Один ( D c ̸ = 0)	3 . 8283	3 . 8220	3 . 8149	3 . 8066	3 . 7985	3 . 7916
Два ( D c ̸ =0)	4 . 0683	4 . 0614	4 . 0537	4 . 0446	4 . 0356	4 . 0281
Один ( D c = 0)	0 . 1516	0 . 1511	0 . 1505	0 . 1499	0 . 1492	0 . 1487
Два ( D c = 0)	0 . 1733	0 . 1727	0 . 1720	0 . 1712	0 . 1704	0 . 1697

Из сравнения рисунков 6 (D c = 52.8 м ² /с) и 10 (D c = 0) следует, что картина трехмерного движения концентрационного фронта качественно ничем не отличается от двух решений одномерной задачи, представленных на рисунке 2. На входе в изделие, как в случае нелинейной, так и классической диффузии, имеем узнаваемый логарифмический профиль зависимости концентрации кремния от радиальной координаты. Однако выпуклость распределения при классической диффузии мала (Рис. 10) и является исключительно следствием схождения линий тока к оси симметрии ввиду топологической особенности задачи. В то же время при рассмотрении нелинейной диффузии, помимо этого механизма, происходит дополнительный конвективный массоперенос паров, вследствие этого профиль распределения на входе в изделие при учете нелинейности процесса оказывается практически в 25 раз круче характерного для классической диффузии.

( a )

0.6

S0.5-

Рис. 10. Поле концентрации в стадии установления в горизонтальном сечении на высоте h = 20 см в случае классической диффузии; сегмент с одним тиглем ( а ), сегмент с двумя тиглями ( б )

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

Работа выполнена в рамках крупного научного проекта при финансовой поддержке Минобрнауки России (соглашение № 075-15-2024-535 от 23.04.2024).