автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Динамическая адаптация во многофронтовых задачах Стефана

На правах

Королева Ольга Николаевна

ДИНАМИЧЕСКАЯ АДАПТАЦИЯ ВО МНОГОФРОНТОВЫХ ЗАДАЧАХ СТЕФАНА

специальность 05. 13. 18 Математическое моделирование, численные методы и комплексы программ.

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Москва 2006

Работа выполнена на кафедре математического и компьютерного моделирования Московского гуманитарного университета

Научный руководитель:

доктор физико-математических наук, профессор В.И. Мажукин

Официальные оппоненты:

доктор физико-математических наук, профессор О.С, Мажорова

доктор физико-математических наук, профессор В.А. Гасилов

Ведущая организация: Институт Общей физики им. А.М. Прохорова

Российской Академии Наук

Защита состоится " " 2006 г. на заседании специализированного совета

К 002.058.01 при Институте Математического Моделирования РАН по адресу: 125047, г. Москва, Миусская пл. д 4-А.

С диссертацией можно ознакомиться в библиотеке ИММ РАН. Автореферат разослан "/X " ОС.Ти? 2006 г.

Ученый секретарь диссертационного совета кандидат физико-математических наук П '¡¡Р&'^Р' Н.Г. Прончева

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы

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

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

Цель работы

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

Научная новизна

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

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

Научная и практическая ценность

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

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

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

Апробация диссертации

Материалы диссертации докладывались на Международной конференции "Laser Processing of Advanced Materials and Laser Microtechnologies", (Москва, 22-27 июня 2002); на Международном научном семинаре "4th International Workshop on Fundamentals of Ablation with Short Pulsed Solid State Lasers 2004", (Hirschegg, Austria), Stuttgart, Germany: Forschungsgesellschaft for Strahlwerkzeuge FGSW, 4-6 February 2004; на Первом, Втором и Третьем Международных научных семинарах "Математические модели и моделирование в лазерно-плазменных процессах" (Москва, 28 января - 1 февраля 2004; 25 - 29 января 2005; 31 января - 4 февраля 2006).

Публикации

Основные результаты по теме диссертации опубликованы в 7 работах, указанных в конце автореферата.

Объем и структура диссертации

Диссертация состоит из введения, четырех глав, заключения и списка литературы из 118 наименований. Объем диссертации составляет 108 стр., включает 29 рисунков, 2 таблицы.

Содержание работы

Во введении обосновывается актуальность темы диссертации, сформулирована цель работы и приведено ее краткое содержание.

Первая глава, ввиду большого объема публикаций по адаптирующимся сеткам, посвящена обзору и анализу литературы за последние 30 лет. Обсуждаются основные проблемы, возникающие при построении, генерации и адаптации расчетных сеток для задач математической физики. Отмечается, что в последнее время область вычислительной математики, связанная с конструированием расчетных сеток, выделилась в самостоятельный раздел науки со своими методологией, понятиями, определениями, методами и классификацией. Выделены основные классы задач, требующие применения адаптивных сеток. Выделены основные направления расчета, генерации и перестройки расчетных сеток. Проводится краткий обзор существующих в настоящее время методов генерации и адаптации расчетных сеток для стационарных и нестационарных задач математической физики. Рассмотрена проблема согласования движения узлов разностной сетки и с этой точки зрения показано место метода динамической адаптации в ряду существующих методов построения адаптирующихся к решению сеток. В заключительной части главы рассматриваются методы решения задач типа Стефана, основной особенностью которых является отсутствие явного выражения для скорости движения фазовых границ. Поэтому большинство задач Стефана являются нелинейными и аналитическое решение их известно только для некоторых случаев. Обсуждаются преимущества и недостатки существующих в настоящее время методов численного решения задач Стефана: методов с явным выделением фазовых границ и методов сквозного счета. Отмечается высокая эффективность метода динамической адаптации в таких задачах.

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

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

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

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

х=хо

(а)

vsl

LIQUID

rj(t)

A1

LIQUID

rsi(t)

Ni

' JV

■T-

LASER

X=*L

Г" r^t) r,Ut)

Cr

SOLID I LIQUID I

LASER

X=Xo

(6)

rj(0

A1

rj(0

Cr

-2,3

X=XL

rs3,(t) r,3v(t)

Ni

Рис.1. Схема лазерного воздействия на многослойную мишень.

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

В основу задачи Стефана положено нелинейное уравнение теплопроводности с соответствующими граничными условиями:

at ox ox x0<x<rkv(t),

, H = C (т)т k = sj n = 1,2,3, (1)

к

t>0.

Начальные и граничные условия:

В начальный момент времени /0 температура всех слоев полагалась равной фоновому значению Т0

1 = 10:ТО0,х) = Т0,

Левая граница х = х0 полагалась теплоизолированной

дх

(2)

На межфазных границах плавления - затвердевания х = выписывались два соотношения: дифференциальное условие Стефана и непрерывность температур

ЛХТ^-МАр^ , (4)

дх дх

(Т!(=Т1=Т(=Тт)й,п = 1,2,3, (5)

На контактных границах х = Г" "*' ставились условия идеального контакта

-1Л+/

х=Г"

Л{т)

дт

дх

Л(т)

дг

дх

, п = 1,2,3,

2пП _ ^я

(б)

(7)

На правой облучаемой границе х = /*"(0> = поверхностное испарение описывалось в приближении кнудсеновского слоя:

Лк{т)д^ = А{ТкУО+рк^

дх

(8)

[Рк^ы = Ри(Чь -«)]'", [р* + рки2ки = Ри+ри(ики -и)2}",

[Т0=аг(М)Тк]п,

\ри=ар{М).рн]\

Рн(ЪГ

(9)

(10) (П)

Рн

ЯП

Рн(Тк) = Рьехр

\Ть Ти

при

п = 1,2,3, к = эЛ М=1, ат(М)= 0.633, ар(М )= 0.326,

М= 0, ат(М)=схр(М) = 1, где Тк- температура конденсированной среды, ат(М), ар(М) - коэффициенты Крута, М - число Маха, К - газовая постоянная, р„, Р„ - плотность и давление насыщенного пара, Л(т) - теплопроводность, Ср(т) - теплоемкость, Н -энтальпия, Ьт,Ьи- теплота перехода, соответственно, при плавлении и испарении, ики, (к = - скорости движения фронтов плавления-затвердевания и испарения, А(Т) - коэффициент поглощения, р(Т) - плотность, Р - давление, и -скорость звука, Ть, Тт, - температуры кипения и плавления.

Для численного решения математической модели (1)-(11) использовался метод динамической адаптации, в основу которого положен переход к произвольной нестационарной системе координат с переменными (д,г),

принадлежащей некоторому расчетному пространству Оя г. Якобианом такого

и, Я V

преобразования является функция

V дх _

— =--. Переход к произвольной

Р дц

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

В результате перехода к произвольной нестационарной системе координат исходное дифференциальное уравнение (1) трансформируется в (12) и дополняется уравнением обратного преобразования (13).

д(у/Н)= д{Н<2) дУУ дт дд дд

д у ~дт'

дд

Чб<Ч<Гко-

а 'дх

к дд р.

т>0, п = 1,2,3,

,н=срт

к

к=5,е.

ж=-

л(т)рдт V дд.

(12)

(13)

Частные производные зависимых переменных выражаются следующим образом:

& д? у дд

ск

у дд'

д __ р д р д дх2 у/ у/ дд

где — = -— - скорость движения нестационарной системы координат. 2 -дт р

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

После определения конкретного вида функции ¡2 уравнение (13) используется для построения адаптирующейся сетки. Его разностный аналог описывает динамику узлов сетки, а функция Q осуществляет контролируемое движение узлов сетки, согласованное с динамикой искомого решения. Согласование достигается функциональной зависимостью функции от искомого

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

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

Учёт в функции преобразования зависимости от температуры Т(1,х) позволяет концентрировать узлы в области больших градиентов решения. Исходя из предположения, что высокая теплопроводность металлов не способствует формированию температурных областей с большими градиентами, можно считать доминирующей особенностью задачи Стефана наличие подвижных фазовых границ. В этом случае можно исключить зависимость Q от Т^.х), а учесть только

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

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

где И - свободный параметр, имеющий смысл коэффициента диффузии.

С учетом (14) уравнение обратного преобразования (13) примет вид ду дО д ^дш

— =--— = —£>—г- (15)

дт дд дд дд

При достаточно большом значении коэффициента В любое возмущение функции у/, возникающее на одной из границ или внутри области, будет равномерно распределяться вдоль пространственной переменной q, что после разностной аппроксимации дифференциальной модели обеспечит квазиравномерное распределение узлов сетки в физическом пространстве Г2Х, на каждый момент времени.

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

В результате получено окончательное выражение для коэффициента Э, связывающее его величину с параметрами исходной задачи - геометрическим размером области ¿(() и скоростями перемещения границ и0,иг:

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

Для численного решения в расчётном пространстве г вводилась расчётная сетка с множеством узлов У, в соответствие с которой, в физическом пространстве

£2Х) вводилась расчётная сетка

Расширенная система дифференциальных уравнений (12) - (13), дополненная начальными и граничными условиями, аппроксимировалась семейством консервативных разностных схем. Полученная система разностных уравнений линеаризовалась, а затем решалась методом раздельных прогонок с итерациями.

С помощью построенной математической модели (12) - (13) исследовалась динамика нагрева и фазовых превращений в трехслойной мишени, состоящей из алюминиевой подложки большой толщины ((л{ - 0.5+1 см), двух тонких покрытий никеля и хрома ({щ = {Ст = 10-50 мкм), нагреваемых лазерным импульсом большой

Рис. 2 а,б. Скорости плавления о"; алюминия, никеля, хрома и испарения и<0 верхнего слоя: а) в трёхслойной мишени А1+№+Сг, б) в трёхслойной мишени А1+Сг+№.

длительности (Т( = 1-5 с) с прямоугольным временным профилем и относительно невысокими значениями интенсивности (О = 3-105 - 3-104 вт/см2). Исследовались две ситуации плавления и испарения трехслойной мишени. В первой, рис. 1(а) в качестве верхнего легирующего слоя использовался хром, а во второй, рис. 1(6) — никель. В обеих схемах (рис. 1) лазерное излучение падает на поверхность верхнего покрытия и частично ею отражается. На рис. 2 показаны скорости плавления алюминия, хрома и никеля и испарения верхнего слоя по обеим схемам (рис.1). Видно, что первым в обоих случаях плавится алюминий и расплавляется полностью до начала плавления никеля, затем плавится никель. Хром, в силу своих тсплофизических характеристик плавится последним. Причем, на схеме где он является облучаемым слоем (рис. 2а), еще до начала его плавления, заметным становится, начиная, примерно с ¡~3с процесс поверхностной сублимации и1м тах ~ 7.6-10~3 см/с, т.е. испарение из твердой фазы. Плавится хром, испарившись на 92,8% и, незадолго до окончания импульса полностью испаряется (2а). Такого пе происходит на схеме 16, когда верхним слоем мишени является никель (рис.2б). Динамику процесса плавления каждого из слоев удобно характеризовать с помощью функций у/!, у/1 показывающих во сколько раз изменяются пространственные размеры твёрдой и жидкой фаз под влиянием движения фронтов плавления. На рис. 3 представлено пространственное распределение функций угг, на несколько моментов времени. В начальный момент / = 0 предполагается, что пространственные размеры физического и расчетного пространств совпадают, поэтому у/;{0,х)= 1. Это соотношение сохраняется вплоть до начала плавления

алюминия, 1 = 1.24-КГ4 с, рис. За, когда область твёрдой фазы начинает убывать у/,(11х)<1 до полного исчезновения = а область жидкой фазы возрастать

до своего максимального значения $/<Л1(?,х)»Р.2, рис. 36. В дальнейшем функция ¥е.ы(1'х)> вплоть до начала кристаллизации, остается неизменной, так как = 0.

Здесь же на рис. 36 показан процесс плавления никеля: (*/т(/,*)« 1, (?,*)»./. Для хрома 1, так как процесс плавления

ещй отсутствует. Распределение ^(/.х), на момент / = 1бс, рис. Зв,

соответствует стадии полного расплава алюминия и никеля и процессу плавления хрома, характеризующегося соотношениями у/1С,(/,х)«1, у/(Сг(/,х)» 1.

Важным представляется определение, в результате моделирования, диапазона интенсивности и длительности лазерных импульсов, при которых происходит полное удаление верхнего слоя.

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

0,5 0,4

v.. V,

100-1 80 60 40 20 0

(6) t = 2,78 с

At liquid

0,2 0,1 0,0 x (см)

Ni liquid j Cr solid

Ni solid i

0,010 0,008 0,006 0,004 0,002 0,000

X(CM)

80 60 40 20

(вИ= 3,6 с

Al liquid Ni liquid 1 ! Cr

! ; liquid

i ! Cr

| J^solld

■---г--■-1— ~~u—

0,010 0,008 0,006 0,004 0,002 0,000

x (cm)

Рис.За-в. Пространственно-временные распределения функций у.

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

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

г = и $ = 7 = »7(х,у).

Это преобразование отображает область произвольной формы Пху физического пространства с координатами (х,у) на параметрический квадрат в плоскости криволинейных координат (¡; ,г|) (Рис. 2). Задачей генерации сеток для двумерных задач является установление связи между точками (£ ,т|) регулярной расчетной области и точками (х, у) области С2ху (Рис. 2).

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

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

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

осуществляется с помощью специально выбранных функций Р\ и Р2 в правой части уравнений Пуассона.

К построенным таким образом криволинейным расчетным сеткам предъявляется ряд специальных требований:

1. Желательно, чтобы расстояния между соседними узлами несильно отличались между собой.

2. Ячейки сетки не должны вырождаться, т.е. углы в элементарной четырёхугольной ячейке не приближались бы к "О" или "я",

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

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

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

концентрации узлов на отрезке воздействия источника энергии.

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

На этапе корректор с помощью численного решения уравнений диффузионного типа производится плавное перераспределение узлов вдоль каждого из направлений.

Частичная адаптация исходной сетки обладает тем достоинством, что позволяет использовать физические особенности задачи при начальном задании фиксированных координатных линий.

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

Математическая формулировка классического варианта нестационарной двумерной задачи Стефана сводится к квазилинейному уравнению теплопроводности

дН _ дЩ дЖ2 | дг дх ду

,к = з,£, (17)

Нк=срРкТ, (Щ)к=-лк(т)^-, {1У2)к=-?.к(Т)^

в двух подобластях ОД/) и П; (/) произвольной области О^,, разделенных заранее неизвестной подвижной границей На межфазной границе Г5Д?)

выполняется дифференциальное условие Стефана и предполагается постоянство температуры

\¥1 = IV/, Т„. (18)

Здесь индексы п и г обозначают нормальную и тангенциальную компоненты, индексы у и I обозначают принадлежность вещества к твердой и жидкой фазам, Тш, — температура и теплота плавления (кристаллизации), и1(- скорость движения границы раздела фаз. На границе дО^ области = йД^иО,^) задаются краевые условия вида

№'й\а = /' (19)

где IV = , 1¥2 ) - вектор теплового потока, Я - внешняя нормаль к сОХ},, /- заданная на 8С1ху функция.

Учет испарения осуществляется в рамках однофазного варианта задачи Стефана и характеризуется наличием подвижной границы раздела фаз Г<ц(7): жидкость-пар в области П = Процесс развитого

поверхностного испарения описывается с помощью трех законов сохранения массы, импульса и энергии на этой границе (8) - (11).

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

Решение совмещенного варианта задачи Стефана состоит в определении температурных полей Те{() и скоростей движения и ики фазовых

фронтов (/), Г*„(/). Алгоритм численного решения реальных задач Стефана, описывающих процессы испарения, плавления - затвердевания, естественным образом распадается на два качественно отличающихся этапа. Первый -предназначен для определения температурных полей в области с одной подвижной границей ГАи(г). Он учитывает поверхностное испарение и охватывает стадии нагрева твердой фазы вплоть до равновесной температуры плавления Тт, с последующим охлаждением до начальной температуры 7о после окончания процесса затвердевания. На втором этапе определяются температурные поля в двух подобластях и связанная с ними скорость движения и!( межфазной

границы Г5Д<).

Введем некоторое расчетное пространство в котором определена

произвольная нестационарная криволинейная система координат (£,7, г). При переходе к нестационарной системе координат

Якобианом такого преобразования является функция J

В результате преобразования к нестационарной системе координат исходное дифференциальное уравнение (17) в новых переменных принимает вид

где

где

Щ =

у[дпд{ а^дгг)' 2 |Д дчд£ + д£дп)'

и дополняется уравнениями обратного преобразования

1

дт р\к'

(22)

(21)

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

Так как решение полной задачи Стефана состоит из двух качественно отличающихся этапов, то для построения расчетных сеток используются два вида отображения. Оба преобразования отображают область произвольной формы Пху физического пространства на прямоугольник Д., в плоскости криволинейных координат (£,77). При этом в расчетном пространстве взаимно ортогональные координатные линии образуют равномерную сетку с прямоугольными ячейками. Решение задачи (17)-(19) и (8)-(11) на этапе нагрева-охлаждения без фазовых превращений ничем не отличается от решения типичных нестационарных задач математической физики. На этом этапе задача решается на неподвижной сетке, которую можно построить до начала расчетов с помощью с помощью генератора сеток, действующего по двухэтапной схеме предиктор - корректор, описанной в третьей главе.

Второй этап решения связан с зарождением и распространением фазовой границы Благодаря использованию нестационарной системы координат в

расчетном пространстве линия фазового перехода Г3Дг) оказывается неподвижной и совмещенной с одной из координатных линий. С помощью граничных условий (18) на этой границе вычисляется поток вещества = -р^",, определяющий переток массы из одной подобласти в другую. В физическом пространстве сетки в обеих подобластях перестраиваются на каждом временном слое в

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

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

д2х

2+А/

д2х дг]2

(23)

+ Д

аУ 'дп2

(24)

где коэффициенты зависят от скорости движения узлов фазового фронта и

степени трансформации областей.

На втором этапе решения задачи Стефана, связанном с зарождением новой подобласти С1( и распространением фазовых границ Г5,(г) иГ4ц(0> построение подвижных расчетных сеток осуществлялось посредством численного решения нестационарных уравнений:

дт *д%2 пдц2 (25)

дт # д$2 * дг)2 (26)

Для определения коэффициентов диффузии используются следующие формулы:

[ле

К

41

у (¿ел

/ им

1

где

М

пи*)'

М

О лцО

_ Л?¡/¡А) -

Здесь Ь, ,А£ ^ .Ь^.М^ - длины координатных линий и длины

ребер ячеек в физическом пространстве в момент времени Г и в начальный момент времени, - длины сторон ячеек в расчетном пространстве. Численное

решение уравнений (25), (26) с функциями (23), (24) позволяет определять координаты узлов сетки на каждый момент времени.

При конечно-разностной аппроксимации преобразованной к нестационарной системе координат исходной задачи (17) - (19), (8) - (11) дополненной уравнениями (25) - (26) в области х [0,/0] пространства обобщенных

координат была построена прямоугольная сетка а> с шагами т,Лт} соответственно:

1 = 0.....1-1,

1к+, = = О.....К -1, т^ = + Ат>|

В узлах сетки определялись функции х]к,у1к>6/д.6г,(*> а в серединах ячеек функции Т/+1/2Л+1/2, ^¡+1/2М1р-Нм/2М1/2- к серединам ребер ячеек были отнесены функции '

Аппроксимировалась исходная дифференциальная задача (17) - (19), (8) -(11) интегро-интерполяционным методом с помощью неявной разностной схемы.

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

На рис 4 приведены результаты моделирования лазерного воздействия с

интенсивностью б = 105 и длительностью т 1 ~ 0,3 сек гауссовского импульса см

на свинцовую мишень треугольной формы. В области Оху была построена сетка

19x18 с максимальным сгущением узлов к источнику. Рис. 4а и 4в соответствуют начальной стадии нагрева с сильным сгущением сетки в зоне воздействия. Неподвижная сетка рассчитана до начала расчетов и используется для определения температурных полей в твердом теле при Ттах < Тт. На втором этапе, который начинается с достижения на поверхности температуры плавления Г„, процессы описываются двухфазной задачей Стефана. Для введения новой (жидкой) фазы допускается перегрев облучаемой поверхности на 0.1К. Из соотношения энергии перегрева и теплоты плавления Ьт определяется толщина жидкой фазы, которая

составила примерно 10~7 м и в области Ох у вводится новая подобласть Пх — и П(. С этого момента в расчетах используется перестраиваемые на

каждом временном слое сетки с общим числом узлов 19x23, из них 19x6 относятся к новой фазе (рис. 4г - 4е). На этих рисунках показаны изменения сетки в разные моменты времени по мере образования канала на стадии плавления. Отслеживается движение фронтов плавления и испарения Г(и в различные моменты времени и соответствующая деформация подобластей пропорционально росту скоростей и!( и и(и. Рисунок 4г отражает момент времени начала формирования канала и возникновения двух фронтов: плавления и испарения. Деформация подобласти П1 усиливается при достижении на поверхности температуры равновесного кипения, при которой испарение становится заметным и скорость испарения и[и увеличивается (рис. 4д - 4е). Финальная стадия формирования канала в жидкой фазе показана на рис. 46 и 4с. На

всех рисунках (4а - 4е) отслеживается адаптация расчетной сетки ко всем особенностям решаемой задачи.

у.«

мша ою

НШ8

Рис. 4 а, е. Расчетная сетка в различные моменты времени.

Выполненное моделирование свидетельствует о высокой эффективности метода динамической адаптации для рассматриваемого класса задач лазерной обработки материалов.

В заключении сформулированы основные результаты работы.

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

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

испарения, происходящие в многослойных мишенях под действием лазерного излучения.

3. Предложен алгоритм численного решения многофронтовых нестационарных двумерных задач Стефана с явным выделением межфазных границ в произвольных областях. Разработанная методика использовалась для моделирования процессов лазерного плавления и испарения в двумерных постановках. Исследованы различные режимы воздействия на мишени сложной геометрии.

Публикации по теме диссертации.

1. V.I. Mazhukin, O.N. Koroleva, М.М, Chuiko. Modeling of formation of deep 2D channels in metal targets via laser irradiation. SPIE - The International Society for Optical Engineering. 2002. V 5121, pp. 87 - 97.

2. O.H. Королева, В.И. Мажукин. Математическое моделирование лазерного плавления и испарения многослойных материалов. - ЖВМиМФ, 2006, т.46 №5, с. 910-924.

3. V.I. Mazhukin, М.М. Demin, O.N. Koroleva. Interaction of laser plasma with target surface. In: Berger, P; Dausinger, F.; Fôhl, C. (Eds.): Proc. 4th International Workshop on Fundamentals of Ablation with Short Pulsed Solid State Lasers 2004 (Hirschegg, Austria). Stuttgart, Germany: Forschungsgesellschaft fur Strahlwerkzeuge FGSW (2004).

4. В.И. Мажукин, O.H. Королёва, К. Февьет. Моделирование неравновесной кинетики ионизации в дуговом разряде. Материалы Первого Международного научного семинара LPPM3, Москва, 2004.

5. Vladimir I. Mazhukin, Olga N. Koroleva. Mathematical Modeling of Melting and Evaporation Multilayer Materials. Материалы Второго Международного научного семинара LPPM3, Москва, 2005.

6. О.Н. Королева, В.И. Мажукин, А.В. Мажукин, François Gentils. Математическая модель неравновесной кинетики электрического пробоя молекулярного газа. Материалы Третьего Международного научного семинара LPPM3, Москва, 2006.

7. О.Н. Королёва, В.И. Мажукин. Математическое моделирование лазерного плавления и испарения многослойных материалов. Материалы Третьего Международного научного семинара LPPM3, Москва, 2006.

Напечатано с готового оригинал-макета

Подписано в печать 13.10.2006. Формат 60x84/16. Печ. л. 1,5. Тираж 100 экз. Изд. № . Заказ Л'а 12а?

Издательство Московского гуманитарного университета, 111395, Москва, ул. Юности, д. 5/1.