Содержание к диссертации
Введение
Глава 1. Разностные схемы для решения уравнений тепловой конвекции ..21
1. Уравнения тепловой конвекции в области с неподвижными границами.
1.1. Уравнения в приближении Обербека-Буссинеска 21
1.2. Уравнения в переменных "функция тока - вихрь скорости" 22
1.3. Безразмерный вид уравнений тепловой конвекции в переменных \//, со 24
1.4. Формулы Тома и Вудса приближенного граничного условия для вихря на стенке 26
2. Уравнения тепловой конвекции в области со свободной границей 30
2.1. Кинематическое и динамические граничные условия на свободной поверхности 30
2.2. Условия для вихря и функции тока на свободной поверхности 33
2.3. Уравнение для касательной составляющей скорости на свободной поверхности 35
3. Метод экспоненциальной подгонки для уравнений тепловой конвекции 37
3.1. Аппроксимация дифференциальных операторов 37
3.2. Модифицированный метод неполной факторизации Булеева и схемы расщепления
Выводы по главе 1 55
Глава 2. Исследование конвекции и переноса тепла в жидкости со свободной покоящейся границей .56
1. Физико-математическая модель.
1.1. Уравнения движения 56
1.2. Постановка краевых условий 58
2. Метод решения 59
3. Анализ результатов 61
Выводы по главе 2 74
Глава 3. Исследование конвекции и переноса тепла в жидкости с криволинейными подвижными границами 75
1. Преобразование областей с криволинейными границами в прямоугольные области. Уравнения движения в новых переменных 75
2. Исследование конвекции и переноса тепла в жидкости с подвижной свободной границей жидкость - газ 86
2.1. Физико-математическая модель 86
2.1.1. Уравнения движения 86
2.1.2. Постановка начальных и краевых условий задачи
2.2. Метод расчета 89
2.3. Результаты расчетов 93
3. Решение двумерной задачи Стефана о фазовом переходе 95
3.1. Физико-математическая модель 95
3.1.1. Уравнения движения 96
3.1.2. Постановка начальных и краевых условий задачи
3.2. Метод расчета 97
3.3. Результаты расчетов 102
Выводы по главе 3 105
Заключение 106
Литература
- Безразмерный вид уравнений тепловой конвекции в переменных \//, со
- Модифицированный метод неполной факторизации Булеева и схемы расщепления
- Постановка краевых условий
- Исследование конвекции и переноса тепла в жидкости с подвижной свободной границей жидкость - газ
Введение к работе
Актуальность. Тема диссертации связана с изучением свободноколнек-тивных течений, вызванных эффектами плавучести (конвекция Грасгофа) и термоканиллярности (конвекция Мараигони). Задачи расчета конвективных течений вязкой несжимаемой жидкости возникают при решении многих прикладных проблем, например:
моделировании динамики ледяного покрова водоемов и водочоков;
моделировании процессов кристаллизации расплавов методами Чох-ральского, Бриджмена, бестигельной зонной плавки (БЗП);
исследовании течений в расплавах при пониженной гравитации;
расчете поля температур в многослойной среде с разной теплопроводностью слоев (например, нефть, грунт, вода);
моделировании процессов глубинной геодинамики (таких как образование рельефа Земли, движение материков, формирование структуры недр).
В связи с возрастающей ролью совершенных кристаллов для современной полупроводниковой электроники, оптических кристаллов для элементов лазерных установок, а также увеличением требований к качеству кристаллов, их размерам и объему их производства, становится важным рациональное управление процессом роста кристаллов. Качество будущего кристалла напрямую связано с процессами гидродинамики и тепломассообмена в процессе роста, поэтому велик интерес к исследованиям этих процессов.
При моделировании роста кристаллов необходимо численно решать уравнения тепловой конвекции. Основной трудностью при решении уравнений, описывающих процесс кристаллизации, является наличие большого параметра в правой части уравнения для вихря (число Релея или число Грасгофа). В результате приходится решать уравнение с малым параметром при старшей производной. При интегрировании таких уравнений возникают трудные вычислительные задачи вследствие образования пограничных слоев или внутренних переходных слоев с большими градиентами решений. Это приводит к снижению порядка сходимости и неравномерной сходимо-
сти на равномерных сетках. Для решения возникающих проблем используют различные подходы: повышают порядок точности разностной схемы; в области быстрого изменения решения используют сетку с мелкими шагами. Однако эффективность этих способов уменьшается с увеличением жесткости уравнения. Другой подход состоит в применении метода экспоненциальной подгонки, или метода интегральных тождеств [1-5]. В данной работе метод экспоненциальной подгонки взят за основу для построения вычислительного алгоритма решения уравнений тепловой конвекции.
С помощью построенного в данной работе алгоритма была исследована термоконвективная система с боковым подогревом и верхней неподвижной свободной границей в условиях невесомости и при увеличении силы тяжести. Подобную задачу для условий невесомости ранее решали Бердников B.C. и Талонов В. А. [6]; Zebib А. и др. [7]; Carpenter В.М. и Homsy G.M. [8]. В этих работах представлена достаточно полная теория пограничного елся Марангони для рассматриваемого течения. В настоящей работе приведено дальнейшее развитие этой теории.
Актуальность решения задач о конвективных течениях в условиях пониженной гравитации возросла в последнее время в связи с развитием ракетно-космической техники, сделавшей доступной длительную невесомость. Особенно эффективным использование условий пониженной гравитации может быть в процессах, сопровождающихся фазовыми переходами. К таким процессам относится, например, получение монокристаллов и различных композиционных материалов.
Цель диссертационной работы заключалась в:
построении вычислительного алгоритма решения стационарных и нестационарных задач гидродинамики в рамках уравнений тепловой конвекции в приближении Обербека-Буссинеска, эффективного для широкого диапазона изменения безразмерных параметров;
моделировании течения жидкости с покоящейся свободной границей в условиях невесомости и с нарастанием силы тяготения;
построении математической модели тепловой конвекции в областях со сложной геометрией и ее апробации на задачах, являющихся составляющими общей проблемы моделирования роста кристаллов методом бести-гелыюй зонной плавки.
Научная новизна.
1) Разработан и реализован оригинальный численный алгоритм для реше
ния стационарных и нестационарных задач гидродинамики в рамках
уравнений тепловой конвекции в приближении Обербека-Буссинеска.
2) Осуществлено развитие теории пограничного слоя Марангоии для тече
ния жидкости в нагреваемой сбоку квадратной полости с неподвижной
свободной границей в условиях невесомости и при нарастании силы тя
готения.
Показано, что в промежутке значений параметра Марангоии Ма 0<Ма<102 гравитационная конвекция становится преобладающей, начиная с параметра Релея Ra=I0 .
Найдены диапазоны изменения безразмерных параметров Марангоии и Прандтля, в которых течение не меняется в условиях отсутствия силы тяжести.
С использованием метода наименьших квадратов получены аппрок-симационные выражения для чисел Нуссельта вдоль горячей стенки, как функции от числа Релея и от числа Марангоии
Построены зависимости интегральных чисел Нуссельта вдоль стенок от чисел Прандля, Марангоии и Релея, и тем самым изучена зависимость теплоотдачи от всех безразмерных параметров задачи.
3) Разработана математическая модель для решения задач тепловой кон
векции в областях со сложной геометрией на основе полученных урав
нений тепловой конвекции в новых переменных. Уравнения получены в
результате преобразования областей с криволинейными границами в
прямоугольные области.
Достоверность полученных . результатов подтверждена многочисленными тестовыми расчетами, показавшими хорошее согласие полученных результатов с численными данными других авторов. В расчетах также осуществлялось варьирование параметров сетки и степени ее сгущения в областях больших градиентоп искомых физических величин.
Практическая ценность работы. Разработанный и реализованный автором вычислительный алгоритм и комплекс программ для решения стационарных и нестационарных задач гидродинамики в рамках уравнений тепловой конвекции в приближении Обербека-Буссинеска может быть использован при численном исследовании различных физических и технологических процессов, сопровождающихся конвективными течениями. Полученные результаты для течения жидкости в нагреваемой сбоку полости с неподвижной свободной границей обеспечивают дальнейшее развитие теории пограничного слоя Марангони и могут быть использованы при исследовании течений в условиях пониженной гравитации и при увеличении силы тяжести. Построенная автором математическая модель для решения задач тепловой конвекции в областях с криволинейными подвижными границами позволяет в рамках единого подхода численно исследовать широкий класс конвективных течений в областях со сложной геометрией. Она может быть использована при моделировании роста кристаллов методом бестигельной зонной плавки, расчете поля температур в многослойной среде с разной теплопроводностью слоев (например, нефть, грунт, вода).
На защиту выносятся: Вычислительный алгоритм расчета свободноконвективных течений несжимаемой жидкости в рамках уравнений тепловой конвекции в приближении Обербека-Буссинеска, основанный на методе экспоненциальной подгонки для аппроксимации дифференциальных уравнений и модифицированном методе неполной факторизации Булеева при решении системы конечно-разностных уравнений. Алгоритм имеет первый поря-
док аппроксимации по времени и второй по пространственным переменным.
Результаты численного моделирования течения жидкости в нагреваемой сбоку полости со свободной неподвижной границей жидкость/газ в условиях невесомости и при увеличении силы тяготения с помощью построенного вычислительного алгоритма.
Математическая модель для численного исследования широкого класса конвективных течений в областях со сложной геометрией, построенная на основе полученных уравнений тепловой конвекции в новых переменных в результате преобразования координат.
Апробация работы. Основные научные результаты диссертации докладывались на V Международной Конференции по Моделированию Приборов и Технологий (ICSDV96) (1996, Обнинск); Международной конференции "Математические модели и численные методы механики сплошных сред" (1996, Новосибирск); Втором Сибирском конгрессе по прикладной и индустриальной математике (INPRIM-96) (1996, Новосибирск); Первой Всероссийской конференции по материаловедению и физико-химическим основам технологий получения легированных монокристаллов кремния ("Кремний-96"), (1996, Москва); XVI Международной школе-семинаре по численным методам механики вязкой жидкости (в составе научного мероприятия "Вычислительные технологии - 98") (1998, Новосибирск); Международной конференции "Математические модели и методы их исследования" (1999, Красноярск); Четвертом Сибирском конгрессе по прикладной и индустриальной математике (JNPRIM-2000) (2000, Новосибирск); XVII школе-семинаре по численным методам механики вязкой жидкости (в составе научного мероприятия "Вычислительные технологии - 2000") (2000, Новосибирск); а также обсуждались на семинарах в Институте Вычислительных Технологий СО РАН и Институте Вычислительной Математики и Математической Геофизики СО РАН.
Публикации. По теме диссертации опубликовано 9 работ, из них 3 статьи, 1 электронная публикация, 4 тезисов конференций и I электронные тезисы конференции.
Личный вклад соискателя заключается в обсуждении постановок задач, разработке адекватных численных алгоритмов и методов решения, создании и тестировании алгоритмов и программ, проведении расчетов, участии в интерпретации результатов численного моделирования.
Структура и объём диссертации. Диссертация состоит из введения, трех глав, заключения, серии численных результатов в виде 4 таблиц и 29 рисунков, и списка цитированной литературы из 115 наименований. Полный объем диссертации составляет 119 страниц.
Безразмерный вид уравнений тепловой конвекции в переменных \//, со
Решение многих задач динамики вязкой жидкости осуществляется в переменных "функция тока - вихрь". Начиная с работ Тома А. [90] такой подход называют двухполевым. Основная трудность этого подхода заключается в том, что математическая постановка задачи не содержит граничного условия для вихря на твердой стенке. Существуют различные методы преодоления этой трудности. Наиболее распространенным является способ вычисления граничного условия для со с помощью различных аппроксимаций условий прилипания, которые формулируются д. і 5у/\ „ г, в виде граничных условии на функцию тока: у/\ = —ч = 0, где О - твердая дп \а стенка, п - направление нормали к ней. Уже первые численные расчеты двухполевым методом, представленные в [90], показали, что способ вычисления граничных условии для вихря существенно влияет на сходимость вычислительного метода. Позднее было замечено, что эффективность неявных схем метода дробных шагов [80] для уравнений Нав.ье-Стокса теряется из-за неустойчивости процедуры вычисления вихря на границе. Поэтому проблеме аппроксимации вихря на твердой стенке уделяется много внимания различными авторами (работы [4], [69], [72], [73], [75], [76], [77], [90]-[98]).
Несмотря на то, что это условие первого порядка точности, оно считается очень надежным и приводит к результатам, хорошо согласующимся с результатами, полученными при помощи форм высших порядков граничного условия для вихря [91]. В работе [76], на основе формулы Тома, Вудс Л. предложил граничное условие для вихря второго порядка точности, которое часто используется при решении краевых задач конвекции и теплообмена: ffli= Tw+p2) (L24) Еще одна часто используемая формула второго порядка точности для граничного условия для вихря была впервые предложена Йенсеном В.Г. [99], а позднее Пирсоном Ц.Е. [95]. В [4], [91] и других работах показано, что решения, полученные с помощью формул Вудса и Пирсона, менее устойчивы, чем при использовании формулы Тома. В случае, когда они дают устойчивые решения, результаты получаются менее точными, чем с использованием формулы Тома (хотя условие Тома - первого порядка точности, а формулы Вудса и Пирсона - второго).
Необходимо сделать замечание о возможной переопределенное граничных условий. Для функции тока ставятся 2 граничных условия:
Если рассматривать только уравнение для функции тока в системе (1.11) - (1.13), то каждое из этих условий будет достаточным для нахождения решения. Очевидно, для уравнения Ду/=со нельзя брать оба условия одновременно, т.к. это делает задачу переопределенной. Но условия \[/=0 недостаточно для того, чтобы определить вихрь о) на стенке. Необходимо также использовать градиентное условие (см. вывод формул Тома и Вудса в [75], [76]). Поэтому, за неимением иного граничного условия для вихря, используется градиентное условие, а условие ц =0 берется для уравнения на функцию тока. Это считается единственно правильным распределением данных условий [91]. Для увеличения эффективности неявных схем применяют различные усовершенствования способа аппроксимации. Например, улучшить точность выполнения условий вязкого прилипания на твердой границе удается с помощью процедуры релаксации в граничных условиях для вихря [4]. При решении стационарных задач релаксация позволяет увеличить скорость сходимости итерационного процесса, а при решении нестационарных - значительно увеличить шаг по времени. В данной работе при использовании формул Тома и Вудса оптимальней параметр релаксации определялся путем численного эксперимента. Другой способ усовершенствования способа аппроксимации вихря на стенке был предложен Пирсоном Ц.Е. [95] и развит в работах Грязнова В.Л. и Полежаева В.И. [69], [77]. В работе [69] предлагается метод, заключающийся в том, чтобы граничное условие для вихря ставить не на границе, а внутри основной области, где вихрь определяется соответствующим уравнением. Уравнение для вихря при этом решается во вспомогательной области, расположенной внутри основной; и связь между вихрем и функцией тока переносится с границы внутрь области. (Подробнее алгоритм [69] можно посмотреть в гл.2, 2.) В [69] показано, что решения стационарных задач, полученные с помощью этого метода и метода применения формулы Тома, практически не различаются. Важным преимуществом метода [69] в сочетании с неявными схемами по сравнению с использованием приближенных граничных условий для вихря (формул Тома, Вудса и т.п.) является вычислительная устойчивость, позволяющая существенно увеличить временной шаг при решении нестационарных задач и уменьшить число итераций при решении стационарных. В данной работе для решения стационарных задач тепловой конвекции используются формулы Тома, Вудса и метод расчета граничных условий [69]. Результаты вычислений с использованием формул Тома и Вудса практически не различаются. Метод [69] позволил улучшить сходимость в 2-4 раза по сравнению с методом использования приближенных граничных условий для вихря.
Модифицированный метод неполной факторизации Булеева и схемы расщепления
В приложениях часто возникает ситуация, когда жидкость граничит с газом. Известно, что динамические коэффициенты вязкости газа на 1-2 порядка меньше соответствующих коэффициентов жидкости. Поэтому можно считать, что при движении газа с умеренными скоростями касательные напряжения, возникающие на границе раздела со стороны газа, пренебрежимо малы. Нормальные напряжения с точностью до знака совпадают с давлением газа [84].
Пусть жидкость заполняет некоторую полость (см. рис.3) и у = J(t,x) -уравнение свободной поверхности. Тогда Таким образом, получаем условие непрерывности касательной компоненты вектора напряжений на свободной поверхности (или условие скорость жидкости в направлении нормали к свободной поверхности. Мы будем искать только устойчивые стационарные решения, поэтому будем считать, что vn = 0. аналог двухполевого метода записи уравнении тепловой конвекции в переменных функция тока, вихрь, записанного для свободной поверхности. Эта идея впервые предложена Овчаровой АС. в работе [35] и может быть эффективно использована для решения уравнений в переменных ці, со со свободной поверхностью жидкость/газ. 3. Метод экспоненциальной подгонки для уравнений тепловой конвекции. 3.1. Аппроксимация дифференциальных операторов.
При рассмотрении уравнений тепловой конвекции (1.17) - (1.19) будем применять метод экспоненциальной подгонки к уравнениям для температуры и вихря. Рассмотрим суть этого метода на примере стационарного уравнения для вихря:
В рассматриваемой прямоугольной области введем прямоугольную сетку с узлами (Xj, ) j), i=l,... ,ЛГ; ./ =/,... ,М Шаги сетки Sxt = xi+1 х„ 8) j - у}-ц - у)-. Таким образом, вся область оказывается разбитой на некоторое число непересекающихся контрольных объемов. Дифференциальное уравнение интегрируется по каждому контрольному объему. Сначала проинтегрируем уравнение (1.51) в пределах от х, до xi+]\ dco
Из этого уравнения можно получить аналогичное уравнение для температуры, полагая Рг=1 и правую часть, равной нулю.
Таким образом, с помощью метода экспоненциальной подгонки мы получили конечно-разностные уравнения для температуры и вихря. Как уже говорилось во введении, основными моментами в этом методе являются совместная аппроксимация конвективных и диффузионных членов (традиционно используется раздельная аппроксимация этих слагаемых: центральные или направленные против потока разности для приближения первой производной и стандартная трехточечная аппроксимация второй производной [52], [53], [100]); а также выполнение законов сохранения для конечных контрольных объемов. Таким образом, даже решение на грубой сетке удовлетворяет точным интегральным балансам.
Уравнение для функции тока не содержит конвективных членов, а содержит только диффузионные. Это уравнение аппроксимируем традиционным способом приближения второй производной. На неравномерной сетке эта аппроксимация выглядит так: 1
Итак, мы представили все уравнения тепловой конвекции в конечно-разностном виде. Каждое из уравнений записано на пятиточечном шаблоне и может быть представлено в виде матричного уравнения где Г- вектор решения, q - вектор правой части, [А] - пятидиагональная матрица коэффициентов. В приложениях может возникнуть ситуация, в которой матрица коэффициентов [А] не удовлетворяет условию диагонального преобладания (например, при численном решении уравнений с малым параметром при старшей производной). Основная сложность решения уравнений, не удовлетворяющих условию диагонального преобладания, связана с тем, что многие эффективные численные методы в этом случае прекращают сходиться или сходятся очень медленно [52]. В этой ситуации традиционно вплоть до настоящего времени используется метод факторизации Булеева [101], который позволяет находить решение и при отсутствии диагонального преобладания. В работе [79] предложена модификация метода Булеева, которая тестируется на стационарном двумерном уравнении теплопроводности и оказывается более эффективной по сравнению с методом Булеева. Суть этого модифицированного метода состоит в следующем [79]. Вместо уравнения (1.57} рассматривается модифицированное матричное уравнение: [A+BJT=q+[BJT. Матрица [В] подбирается так, что разложение [А+В] в произведение нижне- и верхнетреугольиых матриц осуществляется за существенно меньшее количество вычислений, чем непосредственное разложение [А]. Пусть [A ]=[A+B]=[L][UJ, где матрицы [L] и [U] являются соответственно нижне- и верхнетреугольными. Матрица [В] состоит из Каждое из уравнений (1.61) решается с помощью алгоритма прогонки [52]. Для линейных уравнений схема продольно-поперечной прогонки является безусловно устойчивой для любых шагов г [80]. В случае нелинейных уравнений это утверждение, вообще говоря, неверно. При численных расчетах, проводимых по схеме продольно-поперечной прогонки, параметр т подбирался экспериментально.
Постановка краевых условий
Также были построены зависимости интегральных чисел Нуссельта от чисел Прандтля, Мараигони и Релея (соответственно рис. 20, 21, 22). Из зависимости Nu от Рг (рис.20) для различных Ма видно, что чем больше число Марангони для одного и того же числа Прандтля (т.е. чем больше термокапиллярная конвекция), тем сильнее теплоотдача. Для Ма от 0 до 102 коэффициент теплоотдачи не меняется с ростом Рг. При Ма=-103; Ю4 и, начиная с Рг=10, он выходит на некоторый постоянный закон, причем для каждого числа Ма этот закон свой. Кроме того, в промежутке Ма от 103 до 10 наблюдается скачок числа Нуссельта, особенно для случая Рг=1. Причиной этого, по-видимому, является возникновение сложной структуры течения и потеря им устойчивости [33], вследствие чего резко увеличивается теплоотдача. Такой скачок теплоотдачи наблюдается и на рис.21, характеризующем зависимость Nu от Ма для различных Рг. Число Нуссельта увеличивается с ростом числа Релея, характеризующего интенсивность гравитационной конвекции (рис.22). В промежутке Ra от 10 до 10 также происходит скачок Nu из-за возникновения сложной структуры течения. На рис. 21, 22 треугольниками показаны функции, аппроксимирующие зависимость числа Нуссельта вдоль горячей стенки от числа Марангони в промежутке Ма от 103 до 104 и зависимость числа
С помощью построенного вычислительного алгоритма на основе метода экспоненциальной подгонки и модифицированного метода факторизации Булеева было осуществлено дальнейшее развитие теории пограничного слоя Марангони для конвективного течения жидкости в нагреваемой сбоку квадратной полости с верхней неподвижной свободной границей. Были получены следующие результаты.
1. Исследовано влияние гравитационной и термокапиллярной конвекции на поле температур и скоростей в широком диапазоне изменения безразмерных параметров - чисел Марангони: 0 Ма 10 , чисел Релея: 0 Ra 106 и чисел Прандтля: 1 Рг со. Показано, что в промежутке значений параметра Ма 0 Ма 10 гравитационная конвекция становится преобладающей, начиная с Ra lO .
2. Найдены диапазоны изменения безразмерных параметров Марангони и Прандтля, в которых течение не меняется в условиях отсутствия силы тяжести.
3. С использованием метода наименьших квадратов, были получены аппроксимационные выражения для чисел Нуссельта вдоль горячей стенки, как функции от числа Релея l(T Ra l(r и от числа Марангони 103 Ма 104.
4. Построены зависимости интегральных чисел Нуссельта вдоль стенок от чисел Прандля, Марангони и Релея, и тем самым изучена зависимость теплоотдачи от всех безразмерных параметров задачи. Глава 3. Исследование конвекции и переноса тепла в жидкости с криволинейными подвижными границами.
1. Преобразование областей с криволинейными границами в прямоугольные области. Уравнения движения в новых переменных. В данном параграфе приводится выполненное автором преобразование криволинейных координат и записываются уравнения тепловой конвекции в новых переменных. Эскиз к горизонтальному варианту бестигелыюи зонной плавки (вверху) и преобразование координат.
Рассмотрим геометрическую постановку модельной задачи моделирования роста кристаллов методом бсстигсльной зонной плавки (рис.23). Выполним преобразование координат так, чтобы границы раздела фаз gi и g2, а также свободная граница жидкость/газ / стали координатными линиями прямоугольной сетки.
Преобразуем уравнения тепловой конвекции (1.14) - (1.16) при а = л12 к новым переменным у, /л на примере уравнения для вихря. После несложных арифметических преобразований конвективные члены преобразуются к виду:
Аналогично уравнению для вихря преобразуются уравнения для температуры и функции тока. В новых переменных у, /л уравнения тепловой конвекции для жидкой фазы в случае, когда сила тяжести направлена горизонтально справа налево, будут выглядеть следующим образом:
Описанный вариант преобразования координат соответствует горизонтальному варианту метода бестигелыюй зонной плавки. Существует еще, так называемый, вертикальный вариант БЗП [12] (см. рис.1). Здесь сила тяжести направлена вертикально вниз. Автором данной работы было сделано преобразование координат и для такого варианта. Искомое преобразование также является суперпозицией двух преобразований, каждое из которых действует на одну пространственную координату. Первое преобразование выпрямляет свободную границу / второе выпрямляет границы раздела фаз gi и g2 Выпишем формулы общего преобразования областей А, В и С.
Исследование конвекции и переноса тепла в жидкости с подвижной свободной границей жидкость - газ
Рассмотрим двухслойную среду В и С (см. рис.1), где В - область жидкой фазы вещества (расплав), С - область твердой фазы, g2(x) -граница раздела фаз, gi=l, f=l. Пусть в начальный момент времени жидкость находится в состоянии покоя, а начальное положение границы раздела фаз: g2=2. Исследуем численно процесс плавления твердого вещества при условии, что расплав подогревается снизу. Будем при этом рассматривать различные начальные распределения температуры в расплаве.
Такой тип задач относится к задачам Стефана о фазовом переходе. В данном случае будем решать однофронтовую задачу Стефана, которую несложно обобщить на случай многослойной среды и, соответственно, многофронтовой задачи Стефана. Как отмечалось во введении, подобная ситуация может возникнуть при необходимости расчета поля температур в многослойной среде с различной теплопроводностью слоев (например, нефть, грунт, вода) или в задачах, связанных с кристаллизацией и плавлением.
Задача будет решаться методом явного выделения границы раздела фаз и точного удовлетворения условий на ней (см. введение). Идея метода основана на декомпозиции сложной системы (многослойной среды) на более простые элементы (модули) при некотором наборе граничных условий [10]. В каждом модуле численный расчет проводится независимо от других модулей системы, а условия и функциональные зависимости на общих границах соединяют эти модули в исходную систему. Задача решается при следующих предположениях:
Итак, пусть В - область жидкой фазы вещества, С - область твердой фазы, g2 - граница фазового перехода. В области В имеет место конвективное движение жидкости, которое описывается системой уравнений тепловой конвекции в приближении Обербека-Буссинеска (1.14) - (1.16). В подобластях, занятых твердой фазой, выполняется уравнение теплопроводности (3.10). Искомые функции - вихрь скорости, функция тока, температура и граница раздела фаз g2.
Будем считать, что в начальный момент времени жидкость находится в состоянии покоя, начальное положение границы g2: g2 = 2, а в начальном распределении температуры в области В присутствует малое возмущение от линейного распределения. Будем последовательно брать следующие начальные распределения температуры: Т= -0,5 при у=3. Вертикальные границы х=0 и х=1 теплоизолированы. Для скоростей (на границе области жидкой фазы) ставятся условия скольжения: нормальная составляющая скорости на границе равна нулю. В терминах "вихрь, функция тока" граничные условия выглядят так: у/=0, со=0 на границе области В.
Идея метода решения основана на декомпозиции сложной системы (многослойной среды) на более простые элементы (модули) при некотором наборе граничных условий [10]. В каждом модуле численный расчет проводится независимо от других модулей системы, а условия на общих границах соединяют их в исходную систему.
Осуществим преобразование координат так, чтобы криволинейная граница раздела фаз g2 стала координатной линией прямоугольной сетки. Воспользуемся формулами преобразования координат (3.16), (3.17). Уравнения тепловой конвекции для области В в новых переменных , ц имеют вид (3.19) - (3.21), уравнение теплопроводности для твердой фазы -(3.23).
Итерационный процесс решения задачи строится следующим образом. Пусть в момент времени t=t известны поля температуры Т , функции к
Основная сложность задачи заключается в том, что коэффициенты уравнений (3.19) - (3.21), (3.23) зависят от положения границы g2, которая, в свою очередь, зависит от температуры, а значит, и от всех искомых фунщий. Таким образом, наряду с нелинейностью по у/, мы имеем нелинейность по g2. Для разрешения этих нелинейностей используется описанный выше итерационный процесс. За начальное приближение V„ берется значение с нижнего временного слоя.
Как уже отмечалось, каждое из уравнений для Тс, Тв, у/, со решается методом стабилизирующей поправки. Алгоритм решения уравнений этим методом описан в 2 главы 3, при решении задачи со свободной границей жидкость/газ. В результате применения схемы стабилизирующей поправки получается система конечно-разностных уравнений, аппроксимирующая исходные дифференциальные уравнения с первым порядком по времени и вторым - по пространству. Эта система уравнений имеет трехдиагональную структуру и эффективно решается методом прогонки [52].
Было проведено моделирование процесса плавления твердого вещества в двухслойной среде. Результаты расчетов с различными начальными распределениями температуры (формулы (3.33) - (3.35)) представлены на рис. 27, 28, 29 соответственно. Малые возмущения в начальном распределении температуры приводят к различным формам поверхности фазового перехода g2 вследствие разных конвективных потоков в области жидкой фазы и, соответственно, разных конечных распределений температуры. В случае начальных распределений (3.33), (3.34) в области расплава образуется два симметричных противоположно направленных вихря (рис. 27, 28). При этом граница раздела фаз gj в своей центральной части либо выталкивается из расплава (рис.27), либо втягивается в него (рис.28). При начальном распределении температуры (3.35) в области расплава образуется один вихрь (рис.29). Расчеты показывают, что на образование двух: вихрей тратится больше энергии, чем на образование одного вихря. Поэтому в случае двух вихрей (рис. 27, 28) твердое тело проплавляется меньше, чем в случае одного вихря (рис.29).
Результаты расчетов сравнивались с результатами [10], где подобная задача также решена методом явного выделения границы фазового перехода, но с помощью преобразования одной пространственной координаты. Результаты [10], в свою очередь, сравнивались с результатами вычислений [4], [15], основанных на использовании метода сквозного счета соответствующих уравнений с разрывными коэффициентами, и показали хорошее совпадение результатов для искомых функций (температура, вихрь, функция тока). Таким образом, полученные результаты показали достоверность построенной модели. Модель может быть использована в дальнейшем для моделирования роста кристаллов методом бестигельной зонной плавки.
