Содержание к диссертации
Введение
ГЛАВА 1. Обзор литературы. 8
Введение. 8
1. Методы построения адаптивных сеток в стационарных задачах. 10
2. Методы построения адаптивных сеток в нестационарных задачах . 13
3. Методы решения задач Стефана. 15
ГЛАВА 2. Алгоритм решения многофронтовой одномерной задачи стефана методом динамической адаптации . 20
Введение. 20
1. Постановка задачи. 21
1.1. Начальные и граничные условия. 23
2. Теплофизические характеристики и параметры . 25
3. Алгоритм решения. 28
3.1. Нестационарная произвольная система координат 29
3.1.1. Начальные и граничные условия. 30
3.2. Функция преобразованияС). 32
3.3. Выбор коэффициента D. 33
3.4. Разностная аппроксимация дифференциальной модели в 35
расчетном пространстве
4. Результаты моделирования. 38
5. Заключение 47
ГЛАВА 3. Генерация неподвижных сеток в произвольных двумерных областях 49
1. Общие положения 49
2. Распределение узлов сетки на границе dQxy физической области 52
3. Построение расчетной сетки на этапе предиктор 56
4. Построение расчетной сетки на этапе корректор 58
5. Примеры построения сеток в областях сложной формы 59
ГЛАВА 4. Моделирование формирования лазерным излучением глубоких двумерных каналов в металлических мишенях . 70
Введение 70
1. Постановка задачи 71
2. Постановка задачи в произвольной нестационарной системе координат 74
3. Функция преобразования Q. 78
4. Разностная схема. 80
5. Алгоритм решения 82
6. Задача о промерзании 84
7. Вычислительный эксперимент. 87
Заключение. 93
Список литературы.
- Методы построения адаптивных сеток в нестационарных задачах
- Теплофизические характеристики и параметры
- Построение расчетной сетки на этапе предиктор
- Постановка задачи в произвольной нестационарной системе координат
Методы построения адаптивных сеток в нестационарных задачах
Построение адаптивных сеток для решения нестационарных задач связано с дополнительными трудностями. Возникновение дополнительных трудностей связано с тем, что нестационарные задачи описываются уравнениями эволюционного типа и особенности решения этих уравнений, такие как большие градиенты, ударные волны, контактные и фазовые границы, разрывы могут возникать, исчезать и перемещаться во всей области определения решения. Поэтому от качества построения сетки, адаптирующейся к решению, зависит не только точность найденного решения, но и сама возможность качественного решения поставленной задачи. Однако построение сеток, опираясь на опыт решения подобных стационарных задач, не всегда представляется возможным. Ряд требований предъявляемых к расчетным сеткам, адаптирующимся к решению, для нестационарных задач значительно увеличивается по сравнению с требованиями к сеткам для задач стационарных. В качестве основных можно назвать эффективность сетки, то есть возможность получения более точного решения с меньшим количеством узлов, автоматический контроль движения узлов, возможность добавления и удаления узлов в ходе решения задачи. В построении расчетных сеток для нестационарных задач дальнейшее развитие получили методы геометрического, алгебраического и дифференциального подходов.
Проблема построения расчетных сеток геометрическими и алгебраическими методами, учитывающими лишь подвижность границ расчетной области, значительно усложняется. В ходе расчета быстротекущих процессов может происходить существенная деформация границ расчетных областей, которую не всегда можно предвидеть [7]. При пересчете разностной сетки в соответствии с новым положением границ на каждом шаге расчетов растет ошибка интегрирования уравнений, используемых для построения сетки. Это приводит к тому, что движение узлов сетки перестает соответствовать физическому смыслу задачи, например, границы области могут практически остановиться, а узлы сетки за счет итерационного процесса двигаться с недопустимыми скоростями. Поэтому для расчета подвижных сеток в [7] применяется сложная процедура поэтапного построения сетки для каждого последующего по времени шага решения задачи, использующая вычисление сдвигов внутренних узлов сетки по интерполяционным формулам с последующей минимизацией некоторого функционала. В геометрических методах [42] вводятся подвижные сетки, жестко связанные с движением границы или с одной из особенностей решения. С помощью таких сеток в [42, 43] выделена головная ударная волна. Таким образом, построение сеток в рамках геометрического и алгебраического подходов для задач с подвижными границами связано с применением различных способов преобразования координат, направленных на получение нестационарной системы координат, которая отслеживала бы положение перемещающихся границ области. Опыт решения нестационарных задач геометрическими и алгебраическими методами [7, 8, 42-43] показывает, что их применение для построения сеток ограничено задачами, для которых есть возможность построения сетки до начала решения, с не изменяющейся во времени формой области и не слишком сильной деформацией границ, то есть задачами, не требующими перестройки сетки в ходе решения.
Значительно шире класс задач, в которых построение адаптивных сеток осуществляется методами дифференциального подхода. Наиболее эффективными методами при решении нестационарных задач, как и при решении стационарных, оказались методы, в которых преобразование координат зависит от искомого решения (автоматическое преобразование координат). Соблюдение этого принципа обеспечивает при решении уравнений эволюционного типа оптимальное распределение узлов расчетной сетки. Построение адаптирующихся сеток будет иметь свои особенности для дифференциальных уравнений различного типа. Основной проблемой дифференциальных методов конструирования сеток является выбор характеристики решения, используемой для управления перераспределением узлов сетки в расчетной области. Как и в стационарных задачах, для управления движением узлов используется такая характеристика как погрешность аппроксимации на каждом шаге [44-46]. Движение узлов [44] задается силами притяжения-отталкивания, которые определяются отклонением локальной погрешности аппроксимации от среднего значения этой величины в расчетной области. В [45,46] расположение узлов расчетной сетки определяется таким преобразованием координат, которое обеспечивает заданный порядок локальной погрешности. Применение такого подхода осложняется громоздкостью расчетов и необходимостью применения регуляризующих процедур при переходе на новый временной шаг.
Развитием эквидистантных схем для нестационарных задач является подход, при котором в качестве управляющей характеристики используется требование равнораспределения по пространству некоторой функции от зависимой переменной [12, 47-49]. Такой подход в отечественной литературе получил название [5] метода весовых функций. Метод весовых функций обобщается на многомерные задачи [50, 51], однако ячейки сетки могут перекашиваться из-за отсутствия контроля их формы, что является недостатком метода.
Теплофизические характеристики и параметры
Температурные зависимости плотности р(т) и коэффициентов теплоемкости Ср(т) и теплопроводностиЛ(т) алюминия, хрома и никеля, используемые в расчётах, брались из справочных данных [94] - [99] и представлены на рис. 2. Вертикальными линиями на них отмечены равновесные температуры плавления и испарения каждого из материалов. Отметим, что теплофизические характеристики всех материалов претерпевают разрыв при переходе через значение равновесной температуры плавления, кроме того, поведение теплофизических характеристик рассматриваемых материалов сильно различается с ростом температуры нагрева.
Хром и никель относятся к тяжёлым металлам, плотность которых в твёрдой и жидкой фазах в 2 - 3 раза превосходит плотность алюминия, Рьн Ра ; РАЇ рис.2(а). Соответственно скорости распространения фазовых фронтов будут сильно различаться (обратно пропорционально плотности вещества).
Теплоемкости материалов, рис. 2(6), в твердой фазе представлены температурными зависимостями, а в жидкой — константами. Наиболее теплоемким является алюминий. При комнатной температуре его теплоёмкость в 2 раза превосходит теплоёмкость хрома и никеля. С ростом температуры это различие уменьшается. Хром и никель относятся к классу слабо теплопроводящих металлов. Их теплопроводность, рис. 2(в), в 2.5 - 4 раза ниже, чем у алюминия. Обращает на себя внимание наибольший скачек теплопроводности у алюминия (примерно в 3 раза) при переходе фазовой температуры Тт.
Одной из качественных характеристик теплофизических свойств материалов является температуропроводность а = Л/Ср, рис. 3, характеризующая глубину теплового влияния и скорость теплопереноса.
В твёрдой фазе температуропроводность алюминия в 3 - 5 раз превосходит коэффициент температуропроводности хрома и никеля. В
Температуропроводность алюминия, никеля, хрома жидкой фазе это различие достигает более чем 10-ти кратного значения. Таким образом, в процессе лазерного воздействия алюминий будет быстро прогреваться на большую глубину, в то время как у хрома и никеля будут разогреваться до больших температур только приповерхностные слои, что найдёт своё отражение в динамике процессов плавления-затвердевания и испарения материалов многослойной мишени.
В таблице 1. приводятся значения теплофизических параметров, характеризующие свойства рассматриваемых материалов многослойной мишени. Из рассматриваемых материалов алюминий имеет самую низкую температуру плавления Тт=933К, а хром является наиболее тугоплавким Тт = 2133К. Обращает на себя внимание необычно близкое расположение
Для численного решения математической модели (2.1)-(2.11) использовался метод динамической адаптации [59, 63, 93], в основу которого положен переход к произвольной нестационарной системе координат с переменными (q,т), принадлежащей некоторому расчетному пространству OqT. Переход к произвольной нестационарной системе координат приводит к тому, что неизвестными являются не только сеточные функции Т/, но и координаты узлов сетки х{. Для их определения используется уравнение обратного преобразования, представляющее собой дифференциальное уравнение в частных производных. Уравнение обратного преобразования, составляется таким образом, чтобы скорость движения узлов зависела от динамики решения уравнений, описывающих физические процессы.
Переход из физического пространства Qxt в расчётное Oqt осуществляется с помощью замены общего вида х = % (q,r), ( = т, имеющей обратное невырожденное преобразование q = (p (x,t),r = t. Якобианом такого преобразования является функция — = . Частные производные Р дч зависимых переменных выражаются следующим образом [59, 63, 92, 93]: д___д_ Qd_ dt дт у/ dq дх у/ dq д2 = р д р д дх2 у/ dqy/ dq ск О п где — = -— - скорость движения нестационарной системы координат. (J дт р функция преобразования, заранее неизвестна и подлежит определению. В контексте данной проблемы функция Q имеет смысл потока вещества, Q = -ри. В новых переменных (q, т) уравнение (1) запишется в виде:
Построение расчетной сетки на этапе предиктор
На этапе предиктор строятся начальные сетки, отвечающие в большей или меньшей степени трем сформулированным ранее критериям. Примеры начальных сеток для ряда связных областей сложной геометрии приводятся на рис. 13-19. Построенные сетки достаточно хорошо сгущают узлы в определенной области, однако требование равномерности и невырожденности ячеек выполняется не вполне, так как часть ячеек имеет вытянутую форму и острые углы (рис. 13а - 19а), что хорошо видно на выделенных фрагментах начальных сеток. Такие сетки могут не обеспечивать в полной мере точность аппроксимации дифференциальных операторов, входящих в дифференциальные уравнения и хорошую обусловленность матрицы разностных уравнений. Поэтому помимо приведенных выше критериев разумно потребовать распределения узлов внутри физической области Пху квазиравномерно по каждому направлению.
Необходимость сгущения узлов делает сетку неравномерной, и чем больше сгущение, тем более неравномерна сетка. Прямой противоположностью этому критерию является квазиравномерность распределения узлов внутри области /2 , таким образом, для построения сетки должен достигаться компромисс между взаимоисключающими критериями. Достичь компромисса при помощи процедуры (3.12) - (3.19) в полной мере не удается, поэтому в работе генератора сеток необходим второй этап -корректор, обеспечивающий более равномерное распределение узлов внутри области Qrv.
На этапе коррекции численное преобразование координат осуществляется посредством решения эллиптической системы дифференциальных уравнений второго порядка, описывающих диффузионный характер движения узлов сетки: Q,= Q2= где Dc, DT],- коэффициенты, определяющие степень «расталкивания» узлов.
Для решения системы (3.20), использовался метод установления [101], согласно которому система (3.20) преобразуется к виду: В качестве начального условия для системы (3.21), (3.22) использовалось распределение узлов, полученное на стадии предиктор. Остановимся подробнее на результатах работы двухэтапного генератора сеток предиктор - корректор.
С использованием двухэтапного генератора (3.12) - (3.19) и (3.20) -(3.22) построены сетки для ряда связных областей сложной геометрии со сгущением узлов в предполагаемой области воздействия источника энергии. Примеры расчетных сеток, построенные для различных областей показаны на рис. 13-19. Буквами А, В, С, D обозначены вершины криволинейного прямоугольника в физическом пространстве /2, соответствующие вершинам параметрического квадрата / в плоскости криволинейных координат {%,rj). Область выделения энергии, к которой производится сгущение сетки, на всех рисунках определена в виде отрезка AD&Q . На рисунках 13а - 19а показаны начальные сетки, построенные на этапе предиктора (3.12) - (3.19). В областях, имеющих форму треугольника (рис. 13a - 15a) в вершинах, противоположных области сгущения узлов, сетка особенно неравномерна. Области неравномерности сетки выделены в виде фрагментов с увеличением. Особый интерес представляет фигура на рисунке 15а, где узлы концентрируются в вершине треугольника. На увеличенном фрагменте в области сгущения, где для дальнейших расчетов ячейки должны быть мелкими, видно, что они сильно вытянуты и скошены, а их размеры значительно больше ячеек, образованных следующими узлами. Такая сетка не обеспечит необходимой точности аппроксимации дифференциальных операторов. Простое увеличение числа узлов повлечет увеличение времени расчетов. Качественное улучшение сетки без увеличения количества узлов обеспечивает второй этап работы генератора - корректор (3.20) - (3.22). Результат его работы показан на рис. 136 - 156. Система уравнений (3.21-3.22) решалась со значениями D(=\ Dl} =10 (136, 146) и Dn =100 (156). В областях треугольной формы сетка стала равномерной, что показано на выделенных фрагментах. В вершине треугольной области (рис. 156) сетка значительно улучшилась: стала более равномерной, ячейки маленькими, пригодными для дальнейшего решения задачи.
На рисунках 16-19 изображены сетки, построенные в области Пху, имеющей форму криволинейного четырехугольника с различной локализацией сгущения узлов. Начальные сетки, построенные на этапе предиктора (3.2) - (3.19) показаны на рисунках 16а - 19а. Сетки сильно неравномерны в области сгущения (16а), в вершинах (16а - 19а), а на сторонах (18а) ячейки перепутаны и выходят за границу физической области П-ху. Все особенности неравномерных сеток показаны на выделенных и увеличенных фрагментах рисунков 16а- 19а. После работы корректора (3.21-3.22) с параметрами D{=1, Dn = 10 сетки в четырехугольных областях стали равномерными, удалось избавиться от вытянутых ячеек с острыми углами (рис. 166,176,196) и от перепутанных ячеек (рис. 186).
Постановка задачи в произвольной нестационарной системе координат
При решении многомерных краевых задач в областях произвольной формы конечно-разностную аппроксимацию дифференциальных уравнений удобнее проводить в криволинейных системах координат. Тип криволинейной системы координат определяется требованиями и ограничениями, накладываемыми на расчетную сетку решением конкретной задачи. Все ограничения можно условно разделить на две категории. Одна из них обусловлена физическими особенностями рассматриваемых процессов, другая связана с формой и геометрическими параметрами области. Зачастую обе категории ограничений оказываются взаимосвязанными. В случае задач с произвольной формой области наиболее жесткие ограничения связаны с аппроксимацией граничных условий. Наиболее приемлемыми для подобных ситуаций оказываются такие криволинейные системы координат, в которых границы области совпадают с отрезками координатных линий. Сопрягающиеся с границами расчетной области координатные линии избавляют от необходимости дополнительной интерполяции при аппроксимации граничных условий любого типа.
Решение совмещенного варианта задачи Стефана состоит в определении температурных полей 7 (/), T((t) и скоростей движения vsC и ики фазовых фронтов TsC(t), Tku(t). Алгоритм численного решения реальных задач Стефана, описывающих процессы испарения, плавления -затвердевания, естественным образом распадается на два качественно отличающихся этапа. Первый - предназначен для определения температурных полей в области с одной подвижной границей Tku(t). Он учитывает поверхностное испарение и охватывает стадии нагрева твердой фазы вплоть до равновесной температуры плавления Тт, с последующим охлаждением до начальной температуры Г0 после окончания процесса затвердевания. На втором этапе определяются температурные поля в двух подобластях Cis, С1( и связанная с ними скорость движения usf межфазной границы Ts((t).
Каждый из этапов имеет свои особенности. Основная особенность первого этапа состоит в возможности появления больших градиентов вблизи зоны нагрева внешним неподвижным источником. Особенности второго этапа связаны с взаимодействием двух подвижных межфазных границ Г5 (/) и rku(t), которые определяются поведением скоростей и"( и и ки. Эти особенности должны учитываться при построении расчетных сеток. Особенности первого этапа для задачи с произвольной формой области можно учесть, воспользовавшись генератором сеток (3.12) - (3.22), описанным в предыдущей Главе. С учетом особенностей второго этапа ситуация обстоит гораздо сложнее. Помимо проблем, связанных с подвижной границей, дополнительные осложнения возникают с изменением характерных размеров фаз. В момент плавления толщина зарождающейся фазы может составлять величину, равную нескольким атомным слоям, то есть примерно 10А. С течением времени толщина новой фазы может увеличиться на 3-6 порядков. То же самое, только в обратном порядке, может
происходить с исходной фазой. В этих ситуациях может оказаться, что расчетные сетки, используемые в начале расчетов, становятся непригодными через определенное время. Для продолжения расчета сетки необходимо коренным образом перестраивать. С учетом изложенных особенностей, построение эффективных расчетных сеток в задачах типа Стефана возможно только на основе адаптивных подходов, при использовании которых адаптивная сетка будет динамически связана с основным решением.
Указанным требованиям в полной мере соответствует метод динамической адаптации [63, 70, 107, 111], в котором движение любого элемента физического пространства отслеживается с помощью произвольной нестационарной криволинейной системы координат. Произвольность выбора означает, что закон движения системы координат заранее не задается, как, например, в лагранжевых переменных, а определяется в ходе решения. В результате движение узлов адаптирующейся сетки зависит от эволюции численного решения физической задачи. В частности, в задаче Стефана динамика узлов связана со скоростью распространения фазовых фронтов. Обратное преобразование для пространственных переменных в методе динамической адаптации осуществляется с помощью двух нестационарных дифференциальных уравнений в частных производных. Тем самым искомые сеточные функции и координаты узлов сетки определяются совместно из решения единой дифференциальной модели.
Введем некоторое расчетное пространство Q , в котором определена произвольная нестационарная криволинейная система координат ( ,;/,г). Предположим, что на каждый момент времени существует невырожденное взаимно однозначное преобразование = (х,у,(), jj = j(x,y,t), г = /, отображающее физическую область произвольной формы С1ху в прямоугольник Q ;/ в плоскости криволинейных координат (,?/)
