М ЕХАНИКА ЖИДКОСТИ И ГАЗА № 2 • 2014
УДК 532.511:536.252
© 2014 г. О. А. БЕССОНОВ, В. И. ПОЛЕЖАЕВ |
КАРТА РЕЖИМОВ И ПРОСТРАНСТВЕННЫЕ ЭФФЕКТЫ КОНВЕКТИВНЫХ ВЗАИМОДЕЙСТВИЙ В ГИДРОДИНАМИЧЕСКОЙ МОДЕЛИ
МЕТОДА ЧОХРАЛЬСКОГО
Представлена обобщенная диаграмма критических чисел Грасгофа в зависимости от числа Прандтля, показывающая зоны возникновения колебания полей течения и температуры в осе-симметричной и трехмерной постановках модели метода выращивания кристаллов вытягиванием из расплава. Обсуждаются структура термиков при высоких числах Прандтля, особенности трехмерной структуры конвекции в зонах стабилизации и смены опасных мод при разных числах Прандтля. Приведены оценки влияния вращения кристалла на структуру полей течения и температуры.
Ключевые слова: тепловая гравитационная конвекция, рост кристаллов из расплава, гидродинамическая модель метода Чохральского, численное моделирование, конвективная устойчивость, конвективные взаимодействия, смена опасных мод, термики.
Колебания температуры расплава в тигле, являющиеся главной причиной неодно-родностей монокристаллов при их выращивании методом вытягивания из расплава, вызываются конвективными неустойчивостями, в которых основную роль играет тепловая гравитационная конвекция, в большинстве случаев находящаяся под действием вращения. Характер неоднородностей существенно зависит от физических свойств, определяющих число Прандтля Рг расплава. При малых Рг характерна полосчатая неоднородность кристаллов, а при больших неоднородности имеют спиралевидную форму. Режимы этой сложной многопараметрической задачи изучаются вплоть до настоящего времени с использованием технологических экспериментов и с привлечением методов математического моделирования, включая специальные тесты [1—4]. В связи с широким диапазоном материалов, применяемых в современной электронике и оптоэлектронике, физические свойства их расплавов существенно различаются, поэтому важную роль имеет рассмотрение режимов конвективных неустойчивостей тепловой гравитационной конвекции на диаграмме чисел Прандтля Рг и Грасгофа Ог. Дополнительный импульс этим исследованиям вызван изучением возможностей улучшения характеристик кристаллов при их выращивании в условиях микрогравитации и поисками соответствующих наземных альтернатив (см., например, [5, 6]).
К настоящему времени разработаны методы расчета трехмерных течений, а также критических чисел возникновения неустойчивостей в двух- и трехмерных случаях при малых и средних числах Прандтля [6—8]. При этом анализ влияния больших чисел Прандтля удалось выполнить лишь в последних работах [3, 4], исключив из рассмотрения вращение. Однако пространственные эффекты в азимутальной плоскости, представляющие технологический интерес, а также более подробные сценарии движения термиков при больших Рг в [4] не представлены.
В настоящей статье приводятся дополнения к результатам параметрических исследований [4], на основе которых построена обобщенная карта конвективных неустой-чивостей на диаграмме Рг—Ог с пояснением особенностей всех найденных режимов, в том числе пространственных эффектов в азимутальной плоскости. Однако итог этих
работ лишь промежуточен, так как особый технологический интерес представляют режимы с вращением. В частности, интересна интерпретация структур, наблюдавшихся в лабораторных экспериментах [9, 10].
Одна из целей статьи — дать, опираясь на опыт предшествующих работ [11, 12], иллюстрацию того, как работает карта режимов для интерпретации результатов моделирования многопараметрических процессов со сложной пространственной структурой.
1. Математическая модель и численный метод. Установка гидродинамической модели метода Чохральского представляет из себя цилиндрический тигель радиуса Rc, заполненный расплавом до высоты H (см. фиг. 1, [3]). В центре к поверхности расплава примыкает диск радиуса Rx, имитирующий торец выращиваемого кристалла. Открытая поверхность расплава предполагается плоской.
Численное моделирование основано на решении пространственных нестационарных уравнений Навье—Стокса в приближении Буссинеска для несжимаемой вязкой жидкости в переменных скорость V, давление p, температура 9
dV + V • (VV) = -Vp + V2V - Grgö
д t
V-V = 0 (1.1)
^Ö + V- (V0) = — V 20
д t Pr
Здесь Gr = gß LA T/v2 — число Грасгофа, Pr = v/a — число Прандтля, L и AT — масштабы длины и температуры, ß — коэффициент температурного расширения, v — кинематическая вязкость, a — температуропроводность, g — вектор силы тяжести. Масштаб скорости определяется как V0 = v/L, масштаб времени — t0 = L /v. Используется цилиндрическая система координат 0 <ф< 2п, 0 < z < H, 0 < r < Rc.
В соответствии с постановкой задачи в [4] геометрические параметры расчета определены следующим образом: Rc = 1, H = 1, Rx = 0.4. На поверхности расплава выполняются условия скольжения: dVJdz = 0, ôV^/ôz = 0, Vz = 0. Для тигля и кристалла может быть задано вращение с угловыми скоростями, соответственно, Rec и Rex.
Используются следующие граничные условия для температуры: нагрев на боковой цилиндрической границе 9 = 1, охлаждение на границе расплава и кристалла 9 = 0, тепловая изоляция на дне тигля d9/dz = 0. На поверхности расплава применяются два варианта условий: с линейным профилем температуры 9 = (r — Rx)/(Rc — Rx) и с тепловой изоляцией d9/dz = 0.
Для решения уравнений (1.1) используется метод конечных объемов. Применяются разнесенные сетки со сгущениями около дна тигля и у поверхности расплава. Уравнения решаются раздельно, с использованием метода проекций и частично неявной схемы интегрирования по времени. Для решения уравнения Пуассона для давления применен экономичный прямой метод с использованием быстрого преобразования Фурье в азимутальном направлении и разложением по базису собственных функций оператора Лапласа в осевом направлении [13]. Расчетная программа распараллелена для многопроцессорных и многоядерных вычислительных систем с общей памятью в модели OpenMP.
Численный метод разрабатывался в соответствии с концепцией прецизионного вычислительного кода, подразумевающей сбалансированность всех элементов расчетной схемы с точки зрения точности, адекватности и вычислительных затрат. При этом не предполагалось использование схем дискретизации и интегрирования по времени с точностью выше второго порядка, совместного решения систем уравнений и прочих дорогостоящих усложнений, поскольку при полномасштабном совместном примене-
Фиг. 1. Сводная карта режимов неустойчивостей тепловой гравитационной конвекции в координатах Рг—Огс в варианте с теплоизолированной поверхностью расплава: 1 — осесимметричный режим моделирования, 2 — пространственный режим, 3 — позиция смены опасных мод колебаний. Характеристики течения: A — стационарное в обоих режимах моделирования, B — нестационарное и неосесимметричное в пространственном режиме, C — всегда нестационарное, D — нестационарное с сохранением осевой симметрии. Характерные картины течения: а — низкие Рг, преобладание глобального механизма конвекции; б — средние Рг, зона стабилизации; в — высокие Рг, неустойчивость типа Рэлея—Бенара с осесимметричными термиками; г — высокие Рг, неосесимметричные термики
нии они многократно увеличивают вычислительные затраты без адекватного улучшения качества расчета, а при применении по отдельности — не достигают ожидаемого эффекта из-за несогласованности с остальными частями алгоритма.
Метод является прямым и не содержит дорогостоящих итерационных этапов. Благодаря высокой эффективности он может использоваться для детального моделирования ламинарных, переходных и турбулентных пространственных течений на подробных сетках. Результаты сравнения с экспериментом показывают, что используемый метод обеспечивает высокую точность и реалистичность моделирования сложных динамических пространственных течений. Развернутое изложение численного метода представлено в [13].
2. Сводная карта конвективных неустойчивостей расплава. Режимы конвективных неустойчивостей в координатах Рг—Ог, показанные на сводной карте (фиг. 1), определяются зависимостью критического числа Грасгофа от числа Прандтля. На диаграмме представлены две зависимости — для осесимметричного (1) и пространственного (2) режимов моделирования. Заметим, что существенную роль играют и граничные условия на поверхности расплава. Зависимость критического числа Грасгофа Огс от Рг для другого предельного случая, линейного распределения температуры на поверхности расплава, приведена в [4].
Кривые на фиг. 1 разделяют зоны с различными режимами течения. Для точек с координатами (Рг, Ог), расположенных на диаграмме выше соответствующей кривой, течение имеет колебательный характер, а для точек, расположенных ниже, колебания отсутствуют. Нахождение таких зон (областей стабилизации) является целью моделирования в рамках данной постановки задачи, а также одной из целей организации процесса вытягивания монокристалла из расплава. В частности, в условиях микрогравитации, когда число Ог пропорционально микроускорению в полете, реализуется стационарный режим течения, что достигалось в некоторых технологических экспериментах. Это приводило к исключению полосчатой неоднородности монокристаллов в тех случаях, когда числа Грасгофа в расплаве в условиях микрогравитации оказывались, в соответствии с диаграммой на фиг. 1, ниже критических [5, 6]. Заметим, что во многих случаях для этого достаточен уровень микрогравитации g/g0 ~ 10-2—10—3. Это выше, чем уровень микроускорений в экспериментах, где требуется обеспечить высокую степень макронеоднородности монокристалла. В последнем случае необходим диффузионный режим, т.е. полное подавление конвекции, что требует значительно более низких микроускорений вплоть до g/g0 ~ 10—6—10-7 [14, 15]. Для выводов такого характера знание карты режимов оказывается весьма полезным.
Земной альтернативой улучшения однородности монокристаллов является снижение интенсивности конвективного течения в расплаве с помощью некоторого управляющего воздействия, например, динамического, теплового или электромагнитного типа, либо изменением геометрии [16].
На выбор таког
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.