CC BY
100
Вычислительные технологии
Том 8, № 5, 2003
О ПОСТРОЕНИИ ОРТОГОНАЛЬНОЙ РАЗНОСТНОЙ СЕТКИ В КРИВОЛИНЕЙНОМ ЧЕТЫРЕХУГОЛЬНИКЕ
Ю. В. Пивоваров Институт гидродинамики им.. М. А. Лаврентьева СО РАН
Новосибирск, Россия
The paper addresses a problem of construction of the orthogonal difference grid in application to the convection calculation problem in the floating zone. The nonorthogonal grid is constructed using transfinite interpolation method. The orthogonalization of the grid is implemented according to the algorithm, developed by V. A. Veretentsev.
Введение
Для получения монокристаллов кремния большого (более 5 см) радиуса используется метод бестигельной зонной плавки в магнитном поле. В вертикальном цилиндрическом образце с помощью источника высокочастотного электромагнитного поля (индуктора) формируется расплавленная зона, удерживаемая между твердыми частями образца силами поверхностного натяжения и магнитного давления. Эта зона медленно протягивается по всей длине образца снизу вверх. Твердые части образца вращаются в противоположных направлениях. Под действием пондеромоторной, термокапиллярной, центробежной сил и силы плавучести в расплаве возникает нестационарная (из-за дестабилизирующего влияния вращения [1]) конвекция, которая на порядок снижает его перегрев и оказывает существенное влияние на величину и форму плавающей зоны.
В работе [1] в предположении осесимметричности задачи произведен расчет электромагнитного поля, формы границ фазового перехода, свободной границы, границы жидкой пленки, поля температуры в образце и поля скоростей в плавающей зоне с учетом всех перечисленных факторов. При этом при расчете конвекции пондеромоторная сила, сосредоточенная в тонком скин-слое, примыкающем к свободной границе, снесена из уравнения импульса в граничное условие для вихря, что позволило не слишком мельчить сетку вблизи свободной границы. Этот прием в [1] удалось провести ценой отбрасывания конвективных членов в скин-слое. Однако согласно оценкам, сделанным в работе [2], конвективные члены в скин-слое имеют величину того же порядка, что и вязкие.
В работе [3] рассмотрена и решена модельная стационарная (ввиду отсутствия вращения) осесимметричная задача о бестигельной зонной плавке в магнитном поле. Однако магнитное поле индуктора здесь учитывалось только как источник внешнего давления на свободную поверхность и не учитывалось как источник электроконвекции в расплаве.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2003.
Для расчета конвекции в плавающей зоне с учетом пондеромоторной силы, действующей внутри расплава, необходимо построить сетку, сильно сгущающуюся как в окрестности свободной границы, где действует пондеромоторная сила и выделяется джоулево тепло, так и в окрестности остальных границ, где образуются пограничные слои. В настоящей работе построена последовательность ортогональных разностных сеток в области, имеющей типичную форму плавающей зоны. Заметим, что применение ортогональной сетки при решении таких задач не является традиционным: в работе [1] использовался метод конечных элементов, а в [3] расчет производился на неортогональной сетке.
1. Построение неортогональной сетки
Пусть в плоскости гОг дан криволинейный четырехугольник АВСД (рис. 1), причем границы АВ и СД однозначно проектируются на ось Ог и заданы уравнениями
АВ : г = /ц(г), г0! < г < г02,
СД : г = /21 (г), г0! < г < г02, а границы ВС и ДА однозначно проектируются на ось Ог и заданы уравнениями
ВС : г = /12(г), г01 < г < г02,
ДА : г = /22(г), г01 < г < г02-
Пусть в плоскости хО;у дан единичный квадрат К : 0 < х < 1; 0 < у < 1, в котором задана неравномерная сетка 0 = х0 < х1 < ... < х^-1 < х^ = 1; 0 = у0 < у1 < • • • < Ум-1 < Ум = 1- Построим взаимно-однозначное отображение К на АВСД. Будем говорить, что сетка £п, введенная на отрезке [а, Ь], подобна сетке хп, если £п = а(1 — хп) + Ьхп, п = 0, N. Аналогично будем говорить, что сетка пт, введенная на отрезке [с, ¿], подобна сетке ym, если пт = c(1 — ут) + ^ут, ш = 0, М.
Введем на отрезках [г01, г02] и [г01, г02] сетки г{п и г^п, подобные сетке хп, а на отрезках [г01,г02] и [г01 ,г02] — сетки г{т и 2^, подобные сетке ут, и положим согласно методу трансфинитной интерполяции (см. формулы (2.4), (2.5) работы [4])
гпт = /12 (г!т)(1 — хп) + /Ы^т^п
Znm = /ii(r11„)(1 - ym) + /2i(r1„)ym, n = 0, N, m = 0, M.
Тогда если
min(/22(z2m) - f12(z1m)) > 0
m
и
min(/21(r2j - /11(rin)) > °,
n
то множество пар (rnm, znm) устанавливает взаимно-однозначное соответствие между точками x = xn, y = ym квадрата K и точками r = rnm, z = znm криволинейного четырехугольника AB CD, причем точки, лежащие на сторонах K, переходят в точки, лежащие на соответствующих частях границы ABCD.
Построим криволинейный четырехугольник, форма которого типична для плавающей зоны [1]. Зададим фронт кристаллизации (линию AB) уравнением
z = Zc — 1.14(1 — r2/R^), 0 < r < Rc,
фронт плавления (линию CD) — уравнением
z = Zf — 0.27(1 + cos(180r/Rf)), 0 < r < Rf
и ось симметрии (линию BC) уравнением
r = 0. Zc — 1.14 < z < Zf — 0.54,
где Rc = 3.4; Zc = 0; Rf = 0.8884; Zf = 0.9020 — безразмерные геометрические параметры (в качестве единицы измерения выбрана величина Г = 1.5 см).
Форма свободной границы (линии AD) вычисляется методом, описанным в [1], при некоторой заданной функции магнитного давления. Так как AD не проектируется однозначно на ось Oz, повернем построенный криволинейный четырехугольник по часовой стрелке на 20 град, построим в нем сетку, как это описывалось выше, и сделаем обратный поворот. При построении сетки в AB CD нужно задать сетку в квадрате K. Она строилась сгущающейся на границах K таким образом, чтобы максимальные шаги уже ортогонализованной сетки в ABCD (см. след. раздел) в направлении нормали к соответствующим границам были примерно в N/15 (M/15) раз меньше толщины погранслоя в окрестности соответствующей границы (для границы AD эта толщина составляет 0.014 см, для остальных границ она примерно в два раза больше [2]). Это приводит к следующим соотношениям:
h12min = xN — xN-1 = 4.275 ■ 10-3/N,
h11 min = x1 — x0 = 16 ■ h12 min,
h21 min = y1 — y0 = 8.55 ■ 10-3/M,
h22 min = yM — yM-1 = 16 ■ h21 min .
Кроме того, при построении сетки в K предполагалось, что величина £n = (xn+1 — xn)/(xn — xn-1) равна £- > 1 при 1 < n < N/2 — 1, 1 при n = N/2 и 1/£+ < 1 при N/2 + 1 < n <
N — 1, причем величины £-, £ + определяются однозначно, например, при N = 30 £- =
1.353652, £+ = 1.650120. Аналогично Пт = (ym+1 — Ут)/(Ут — Ут-1) = П- > 1 при 1 <
m < M/2 — 1, 1 при m = M/2 и 1/п+ < 1 при M/2 + 1 < m < M — 1. При M =
30 п- = 1.554128, п+ = 1.274906. На рис. 2 показана сетка, построенная по описанному выше алгоритму, при N = M = 30.
В
Рис. 2.
2. Ортогонализация разностной сетки
Для построения ортогональной сетки в области Д0, тождественной криволинейному четырехугольнику АВСД используем алгоритм, основанный на минимизации функционала
д К!)’+О' *1 (£ )2+■ (дУ
О0
описанный в [5]. Здесь х(г, г), у (г, г) — функции, отображающие Д0 на К; / — числовая переменная.
Уравнениями Эйлера для функций х(г, г), у(г, г), выражающими необходимое условие существования экстремума функционала, являются уравнения Лапласа, которые после обращения роли зависимых и независимых переменных и деления на (((дг/ду)2 + (дг/ду)2)((дг/дж)2 + (дг/дх)2))1/2 принимают вид
д 2г 5 2г 5 2г д2г д2г д2г ...
— = 0 — = 0 (1) дх2 дхду ду2 дх2 дхду ду2
где
(dr/dy)2 + (dz/dy)2 1
а ^ (dr/dx)2 + (dz/dx)2’ С а’
b (dr/dx)(dr/dy) + (dz/dx)(dz/dy)
^((dr/dy)2 + (dz/dy)2)((dr/dx)2 + (dz/dx)2)
Условие трансверсальности на границе области Do после обращения роли зависимых и независимых переменных принимает вид
dr dr dz dz dx dy + dx dy ’
(r(0’ y)’ z(0’ y)) £ BC’ (r(1,y),z(1,y)) £ DA’
(r(x, 0),z(x, 0)) £ AB, (r(x, 1),z(x, 1)) £ CD. (3)
Оно выражает ортогональность линий x = const и y = const на границе области D0.
На решении задачи а = /, b = 0, так что уравнения (1) после замены x = x, y = /y переходят в уравнения Лапласа. Функции f(x,y) = r(x,y//) и z(x,y) = z(x,y//) осуществляют конформное отображение прямоугольника K : 0 < x < 1, 0 < y < / на D0.
Опишем итерационный процесс построения ортогональной разностной сетки в Д0. Сначала при фиксированном расположении узлов сетки на границе определяются положения внутренних узлов путем приближенного решения разностных аналогов уравнений (1). Это приближенное решение рассчитывается попеременно-треугольным методом при фиксированных функциях а, Ь, с [6].
Далее при фиксированных внутренних узлах ищутся новые положения граничных узлов. Для этого аппроксимируем (2). Например, для границы ВС можем написать
г1т — г0т + (х1 — х0)с0т / Г0то+1 — г0т — г0т — Г0т-^ , х
х1 — х0 (Ут+1 — Ут-1) \ Ут+1 — Ут Ут — Ут-1
х (Г0т+1 — Г0т-Р + (г1т — г0т + (х1 — Х0)сС)т х
(Ут+1 — Ут-1) \ х1 — х0 (Ут+1 — Ут-1)
( г0-1 _ г0 _ г0-1 )) ^г0-1 _ г0-1 ^
X I г0т+1 г0т г0т г0ш-П | уг0т+1 г0т- 1/ __ д (4)
\ Ут+1 Ут Ут Ут-1 / / (Ут+1 Ут-1)
где к — номер итерации. Это равенство имеет второй порядок аппроксимации при условии Ут+1 _ 2Ут + Ут-1 _ О ((Ут+1 _ Ут )2) и условии достаточной близости функций Г„т, ¿„т к решению задачи. Чтобы найти г0т, г0т, нужно определить пересечение прямой (4) с границей ВС. Каждая из четырех границ Д0 задается уравнениями вида г _ гДз), г _ гДз), где гДз), гДз) — параболические сплайны переменной 5. Так что для нахождения очередного приближения положений точек на границе нужно последовательно решать некоторые квадратные или линейные (как в случае границы ВС) уравнения относительно 5, пока не будет найден корень, принадлежащий данному промежутку.
Итак, ортогональная разностная сетка строится в результате многократного чередования описанных итерационных процессов. Заметим, что для ее построения не требуется однозначной проектируемости участков границы на соответствующие оси.
На рис. 3 показана ортогональная сетка в Д0 размерности 120 х 120. Она была построена следующим образом. Сначала была ортогонализована сетка размерности 30 х 30, построенная в предыдущем разделе. Путем интерполяции эта сетка преобразована к размерности 30 х 60 и вновь ортогонализована. Процесс удвоения числа узлов сетки с последующей ее ортогонализацией продолжался, пока не была построена сетка максимальной размерности. Число итераций при решении уравнения (2) с учетом (3) было равно 10. Обозначим через Р число итераций при решении уравнений (1) попеременно-треугольным методом; Q
— число чередований итерационных процессов построения внутренних и граничных узлов
сетки;
біп а
(дг/дж)(дг/ду) — (дг/ду)(дг/дж)
а/ ((дг/ду)2 + (дг/ду)2)((дг/дж)2 + (дг/дж)2) синус угла между координатными линиями;
(/дг/дж — дг/ду) ■ 100% (дг/ду + /дг/дж) ■ 100%
¿1 = ---у / , ч ^ / , ч Г-, ¿2
а/ (/дг/дж)2 + (дг/ду)2
а/ (/дг/дж)2 + (дг/ду)2
— относительные погрешности в выполнении условий Коши — Римана, где / определяется один раз после завершения итерационного процесса построения сетки размерности 120 х 120. Введем также интегральные характеристики:
¿а = У У БіП а^ж^у, ^1 =
к \
¿2^ж^у,
к
^2
\
5|^ж^у, аь
к
\
62^ж^у ■ 100 %.
к
В таблице приведены значения параметров Р, ф, Ж, М, 50,, а2, а при построении
ортогональной сетки.
Здесь первая строка соответствует сетке, построенной в предыдущем пункте. Заметим, что большие значения параметра ф на начальном этапе построения сетки размерности 30 х 30 необходимы ввиду очень малого шага сетки на границе К по нормали к ней. Суммарное время расчета сетки размерности 30 х 30 на ПЭВМ РеП;шт-3 с частотой 1.2 ГГц составляет 8.5 минут. (При необходимости его можно существенно уменьшить, если на начальном этапе временно ввести равномерную сетку в К.) При Р > 128 время счета можно приближенно определить по формуле ¿с = 5 • 10-6РфЖМ с. Из таблицы видно, что
р Я N х М 5« оі 02 оь
0 0 30 х 30 0.88213 50.3 36.6 -
4 4000 30 х 30 0.98945 10.2 11.7 14.4
8 2000 30 х 30 0.99915 3.3 4.0 4.1
16 1000 30 х 30 0.99968 2.4 2.6 2.5
32 500 30 х 30 0.99977 2.2 2.4 2.1
64 250 30 х 30 0.99980 2.2 2.3 2.0
128 125 30 х 30 0.99981 2.2 2.3 1.9
128 125 30 х 60 0.99989 1.8 2.0 1.4
128 60 60 х 60 0.999971 1.4 1.0 0.7
128 60 60 х 120 0.999978 1.2 1.0 0.6
128 60 120 х 120 0.999976 1.1 0.9 0.7
256 60 120 х 120 0.999977 1.1 0.9 0.7
512 1 120 х 120 0.999938 1.1 1.0 0
512 1 60 х 120 0.999946 1.1 1.1 0
512 1 60 х 60 0.999941 1.1 1.1 0
512 1 30 х 60 0.999905 1.4 1.8 0
512 1 30 х 30 0.99988 1.7 2.1 0
при удвоении числа узлов (в каждом из направлений O'x, O'y ) сетки размерности 30 х 30 она значительно приближается к ортогональной (параметры Sa и аь), а при удвоении числа узлов сетки размерности 60 х 60 этого уже не происходит. Вероятной причиной этого является влияние ошибок округления.
Заключительный этап построения ортогональной сетки в D0 — определение l по формуле
1 = // а<Ыу.
K
решение уравнений (1) при a = l, b = 0. c = 1/l с граничными условиями 1-го рода (см. последние пять строк таблицы) и растяжение переменной y в l раз. Значение l = 1.08732 определено по сетке размерности 120 х 120 и затем приписано решению на более редких сетках.
Сравнение результатов расчетов при различных N и M показало, что решение имеет 1.5-й порядок сходимости, хотя все уравнения и граничные условия аппроксимированы с вторым порядком. Это объясняется наличием степенной особенности в тупом угле A области D0. Анализ поведения в окрестности точки A локальных характеристик, соответствующих 12-й строке таблицы, показывает следующее. Параметр sin a терпит минимум, равный 0.9984, а параметр b' = b • 100% — максимум модуля, равный 5.7%, при n =119, m = 1 (x = 0.99996, y = 7,1 • 10-5). (Все характеристики определяются только внутри D0.) При n > 115, m < 2 (x > 0.99977, y < 1.5 • 10-4) имеем sin a < 0.9997, |b'| > 2.5%, а вне этого прямоугольника sin a > 0.9997, |b'| < 2.5%. То есть параметры sin a, b', отвечающие за ортогональность сетки, несколько ухудшаются лишь в малой окрестности точки A. Ухудшение параметров Si и S2, отвечающих за равенство коэффициентов Ламэ в направлениях x, у, более значительно и сосредоточено на бо'льших площадях: max |S1| = 22.6% достигается при n = 118, m = 1 (x = 0.99992. у = 7,1 • 10-5); |S1| в основном больше 6% при n > 76, m < 28 (x > 0.9441. у < 0.01177); max |S2| = 9,0% достигается при n = 119, m = 32 (x = 0.99996, y = 0.018342); |S2| частично больше 6% при n > 92, m < 40 (x > 0.9921. y < 0.04371).
В дальнейшем на построенных сетках предполагается делать расчеты конвекции в плавающей зоне.
Список литературы
[1] MuHLBAuER A., MuizNIEKs A. ET AL. Interface shape, heat transfer and fluid flow in the floating-zone growth of large silicon crystals with the needle-eye technique // J. Crystal Growth. 1995. Vol. 151. P. 66-79.
[2] Пивоваров Ю.В. Параметрический анализ задачи о бестигельной зонной плавке в магнитном поле // Динамика сплошной среды: Сб. науч. тр. / РАН. Сиб. отд-ние. Ин-т гидродинамики. 2000. Вып. 116. С. 142-147.
[3] ОвчАровА А.С. Влияние эффектов плавучести и термокапиллярности на форму свободной поверхности расплава // Вычисл. технологии. 2002. Т. 7. С. 323-328.
[4] ЛисЕйкин В.Д. Метод алгебраической адаптации // Журн. вычисл. математики и мат. физики. 1998. Т. 38, № 10. С. 1692-1709.
[5] Веретенцев В.А. Построение разностной сетки в области с криволинейными границами с помощью конформного отображения // Актуальные вопросы прикл. математики. 1989. М.: Изд-во Моск. ун-та. С. 88-93.
[6] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука. 1978. 592 с.
Поступила в редакцию 10 февраля 2003 г., в переработанном виде — 29 апреля 2003 г.
