Алгоритм численного решения задачи Стефана и его применение к расчетам температуры вольфрама при импульсном воздействии

Полный текст

ОГЛАВЛЕНИЕ 1. Введение 442 2. Постановка задачи 444 3. Численное моделирование 445 4. Результаты численных расчетов 449 5. Выводы 450 Список литературы 451 1. ВВЕДЕНИЕ Явления плавления и затвердевания имеют место в разнообразных природных и технологических процессах, от таяния айсбергов в полярных морях или застывания вулканической лавы до производства мороженого или стали. Главной особенностью этих процессов являются априорно неизвестные, как правило, движущиеся границы раздела фаз. Такие границы принято называть свободными. Для того, чтобы произошел фазовый переход материала из твердого состояния в жидкое, необходимо подвести к образцу тепловую энергию дляразрыва связей, удерживающих ионы в кристаллической решетке. При обратном переходе из жидкого агрегатного состояния в твердое энергия берется из жидкой фазы и направляется на замедление движения ионов и дальнейшей их организации в стабильную решетчатую структуру. Работа выполнена при финансовой поддержке гранта РФФИ, проект № 19-01-00422. © РОССИЙСКИЙ УНИВЕРСИТЕТ ДРУЖБЫ НАРОДОВ, 2021 Эта работа доступна по лицензии Creative Commons 4.0 International https://creativecommons.org/licenses/by-nc-nd/4.0/deed.ru 442 АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ СТЕФАНА 443 Этот интуитивно простой физический процесс описывается математически с помощью классической задачи о фазовых переходах, так называемой задачи Стефана. Решение такой задачи заключается в нахождении поля температур и определении координат границ раздела фаз. Первые подробные исследования задачи Стефана были проделаны в работах французских ученых Г. Ламе и Б. П. Клапейрона [24] и словенского физика и металлурга Й. Стефана [28], имя которого и получила задача. Обе публикации были посвящены процессу эволюции льда в полярных морях. Качественный скачок в изучении задачи Стефана произошел в середине 70-х годов ХХ века, когда задачи со свободными границами выделились в отдельный раздел математики (см. [9] и приведенную там библиографию). В 1973 году Г. Дюво [15] осуществил редукцию многомерной однофазной задачи Стефана к параболическому вариационному неравенству. Это позволило значительно упростить доказательство теоремы о существовании и единственности обобщенных решений, а также получить ряд аналитических результатов о качественных свойствах решений и свободных границ [11-13, 18]. Вариационная формулировка также была установлена и для двухфазной задачи Стефана [17], но полученные результаты по свойствам свободных границ в двухфазном случае оказались слабее, чем в однофазном. Отметим также, что псевдопараболические вариационные неравенства облегчили создание численных методов для решения задач Стефана в многомерных случаях [22, 23]. Помимо аналитических методов, при решении однофазной и двухфазной задач Стефана широко используются методы конечных элементов [4], конечных разностей [30] и интегральные методы [1]. Возможно использование конечно-разностных методов на адаптивной quadtree сетке, что приводит к значительному повышению точности вблизи границ раздела фаз [26]. Некоторые из указанных численных методов входят в специальные математические пакеты. Вычислительные комплексы ANSYS FLUENT, STAR-CD, PHOENICS, развивающиеся в течении последних 30 лет, основаны на эффективных численных алгоритмах, обладают удобным интерфейсом и мощными графическими средствами длявизуализации результатов. В случае существенного вклада не только диффузионного, но и конвекционного переноса тепла определяется поле скоростей в жидкой фазе, получаемое из решения уравнения Навье-Стокса (см., например, [14, 19, 21]). Сложность решаемых инженерных задач ограничена практически только характеристиками компьютера пользователя. Тем не менее, новые поисковые исследовательские задачи требуют построения алгоритмов и программ для моделирования процессов с новыми наборами коэффициентов и их соотношений. В этом случае используютсядва подхода к численному решению. Первый подход состоит в явном выделении расположения границ между фазами. Точное определение ячейки или узла, через которые проходит межфазная граница, часто реализуется с использованием неравномерных и/или динамических сеток [2]. Наиболее известны метод ловли границы в узел пространственной сетки (variable time stepping) и метод выпрямления фронтов (front-fixing method). Второй подход состоит в использовании методов сквозного счета, при использовании которых начальнокраеваязадача решаетсяво всей расчетной области без выделенияобласти фазового перехода. При решении ряда многомерных задач Стефана исключение детального вычисления координат границ между фазами становитсяосновным преимуществом. Реализация таких методов состоит в изменении формы записи исходного уравнения и сглаживания разрывных коэффициентов [3, 5, 6]. Такой подход предпочтителен в случае решения задач, не требующих высокой точности вычисления координат свободных границ. Дополнительных усилий требует минимизация влияния значений параметров сглаживания на погрешность численного решения. В популярных инженерных пакетах программ, основанных на конечно-объемных методах, часто используются методы сквозного счета функций уровня(level set method)и фазового поля (phase field method). Выполнение закона сохранения тепла гарантируется использованием консервативных разностных схем. В настоящей работе представлены метод и результаты моделирования с помощью двухфазной задачи Стефана процесса плавления и испарения вольфрама при облучении его импульсным электронным пучком. Сравнение экспериментально измеренной на установке ВЕТА (см. [29]) зависимости радиуса расплавленной области от времени с расчетными данными, полученными одним из методов сквозного счета, показало хорошее соответствие. Представленная модель позволяет правильно интерпретировать процессы нагрева, испарения и остывания при импульсном тепловом воздействии, что дает необходимые и важные результаты для развития ITER и других экспериментальных термоядерных реакторов. 444 Д. Е. АПУШКИНСКАЯ, Г. Г. ЛАЗАРЕВА РИС. 1. Схема вольфрамовой мишени с вырезанной четвертью. Голубыми стрелками обозначено направление воздействия импульсного излучения лазера. Область расплава вольфрама выделена розовым цветом, Γ - свободная граница между жидкой и твердой фазами, γ - область теплового воздействия. 2. ПОСТАНОВКА ЗАДАЧИ В качестве материала для конструирования диверторных пластин установки ITER предлагается вольфрам, который имеет температуру плавления Tm = 3695 K. На экспериментальном стенде Beam of Electrons for materials Test Applications (BETA), созданного в ИЯФ СО РАН [29], проводилось натурное моделирование нагрева вольфрамовой мишени мощным осесимметричным субмиллисекундным пучком электронов [27]. Образец (прямоугольный параллелепипед) имел размеры 25 мм × 25 мм и типичную толщину 4 мм. В экспериментах происходил мгновенный нагрев пластинки до 9000 K за счет воздействия на её поверхность теплового потока и последующее охлаждение. Распределение плотности мощности нагрева измерялось с помощью рентгеновской визуализации [7]. Именно сочетание высоких плотности мощности воздействия и скорости нагрева привело к необходимости создания новых моделей и программ. Длярасчета распространениятемпературы с нагреваемой поверхности вглубь образца решается задача Стефана в аксиально-симметричной постановке. Поскольку импульсное воздействие лазера на образец продолжалось несколько десятков микросекунд, то образец нагревался на глубину несколько сотен микрон. Поэтому в качестве области моделирования мы выбираем поперечное сечение образца (см. заштрихованный прямоугольник на рис. 1), т. е. рассматриваем двумерную область D := {(r, z) : r ∈ [r0, rmax], z ∈ [z0, zmax]}, (2.1) где r0 = z0 = 0, rmax = 12 мм, zmax = 3 мм. Обозначим через γ часть внешней границы ∂D области D, котораяподвергается тепловому воздействию. Свободную границу между расплавом и твердым вольфрамом обозначим через Γ (см. рис. 1). Тогда задача Стефана принимает вид: ⎧ ∂T ⎪ ⎪c(T )ρ(T ) ∂t ⎪ 1 ∂ = r ∂r ( rλ(T ) ∂T + ∂r ∂ ( λ(T ) ∂z ∂T ∂z в D, ⎨ ⎪ (n, ⎪ 1 T ) ∇ 1γ - W (t, r) N (t, r) = , λ(T ) (2.2) ⎪⎪ (n, ∇T ) = 0 на ∂D\ γ, ⎪ ⎩ T = T0 при t = 0, где T (t, r, z) - температура, c(T ) - удельнаятеплоемкость, ρ(T ) - плотность, λ(T ) - теплопроводность, W (t, r) - мощность теплового потока на γ, N (t, r) - мощность испарения на γ, n - вектор внешней нормали к ∂D, а T0 - начальная температура. АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ СТЕФАНА 445 Распределение мощности по поверхности теплового потока задается формулой W (t, r) = Wmax(t) exp{-Ar2}, A = 0,03088523 мм-2, (2.3) а значение Wmax определено экспериментально. Так как распределение (2.3) близко к нормальному, то вблизи оси симметрии особенностей решениясистемы (2.2) не возникает, как и в других задачах такого типа. На границе контакта твердой и жидкой сред температура предполагается непрерывной. Кроме того, поскольку фазовый переход сопровождаетсявыделением/поглощением определенного количества тепла, то тепловой поток на границе фазового перехода будет разрывен и равен произведению энтальпии фазового перехода на нормальную компоненту скорости движенияграницы раздела фаз. Таким образом, условия на свободной границе Γ имеют вид [T ]1 1 ∂T 1 = 0, λ(T ) = L v , (2.4) 1Γ 1Γ ∂t 1 m n где Lm - энтальпия фазового перехода, а vn - нормальнаясоставляющая скорости движения границы фазового перехода. 3. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ Напомним, что нагреваемый образец имел размеры 25 мм × 25 мм и толщину 4 мм, а область моделированияпредставляла собой поперечное сечение образца размером 12 мм × 3 мм (см. рис. 1). Что касается времени, то численное моделирование продолжалось до того момента, когда было произведено последнее измерение температуры поверхности. Для удобства вычислений, в уравнениях (2.2)-(2.4) был выполнен переход к безразмерным переменным, определенным следующим образом: r∗ = r , λ∗ = r0 λ , ρ∗ = λ0 ρ , c∗ = ρ0 c t , t∗ = c0 t0 = λ0t 0 ρ0c0r2 , T ∗ = T , W ∗ = T0 λ0T0W , r0 где r0, λ0, ρ0, c0 и T0 - это характерные значения параметров, которые приведены в таблице 1. ТАБЛИЦА 1 Параметр Типичное значение Ед. измерения r0 10-1 мм t0 102 мкс λ0 10-1 Вт/мм · K ρ0 10-5 кг/мм3 c0 108 Вт · мкс/кг · K T0 103 K W0 103 Вт/мм2 Точнее говоря, мы, не меняя обозначений, считаем величины в уравнениях (2.2)-(2.4) безразмерными, а для получения результатов в физических величинах после проведения необходимых вычислений домножаем результаты на соответствующие характерные значения параметров. Дляпостроения эффективного вычислительного алгоритма чрезвычайно важное значение имеет тот факт, что задача Стефана допускает обобщенную формулировку (см. [5]), при которой условия (2.4) на свободной границе включаются в первое из уравнений (2.2). Затем коэффициенты полученного уравнения сглаживаются (см. [3]). В результате происходит исключение вычисления координат свободных границ и уравнение теплопроводности во всей области формулируется в единообразной форме с применением дельта-функции Дирака. Точнее говоря, в первом из уравнений (2.2) в множителе при производной по времени появляется дополнительное слагаемое, содержащее теплоту плавления Lm = 51,1 · 105 Вт · мкс мм3 , 446 Д. Е. АПУШКИНСКАЯ, Г. Г. ЛАЗАРЕВА действующую на интервале сглаживания, и уравнение теплопроводности принимает вид: ∂T 1 ∂ ( ∂T ∂ ( ∂T (c(T )ρ(T )+ Lmδ(T, ε)) = ∂t r ∂r rλ(T ) ∂r + λ(T ) ∂z ∂z . (3.1) Здесь множитель δ(T, ε) действует на узком интервале сглаживания [Tm - ε, Tm + ε], ε = 5 K, и вычисляется следующим образом: ⎧ 1 ⎨ , |T - Tm| � ε, δ(T, ε) = 2ε ⎩0, |T - Tm| > ε. m РИС. 2. Графики зависимости от температуры плотности испарения (а), теплопроводности (б), удельной теплоемкости (в), мощности испарения (г), которые использовались при численном моделировании. Заметим, что измерение теплофизических характеристик тугоплавких металлов является сложной задачей. Многие справочники и статьи дают приблизительные или теоретически предсказанные зависимости с оценкой точности менее 10%. Поэтому плотность ρ(T ) (см. рис. 2а), теплопроводность λ(T ) (см. рис. 2б), удельнаятеплоемкость c(T ) (см. рис. 2в) и мощность поверхностного испарения N (t) (см. рис. 2г), стоящие в коэффициентах уравнения (3.1), учитывались согласно зависимостям от температуры материала в диапазоне 300 K � T � 8000 K. Все эти функции имеют разрывы или теряют гладкость при температуре плавления вольфрама Tm = 3695 K. Зависимости от температуры для теплопроводности λ(T ) и теплоемкости c(T ) твердого вольфрама взяты из [16]. Оценки теплопроводности жидкого вольфрама взяты из работ [20, 27]. Зависимость теплопроводности от температуры была скорректирована в процессе моделирования [25]. АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ СТЕФАНА 447 В двумерной области D (см. определение (2.1)) была введена равномернаяпрямоугольная сетка. Размер расчетной области соответствует размерам образцов, но постановка задачи в цилиндрической геометрии предполагает рассмотрение пластинки в форме шайбы. Такое предположение не вносит существенного влияния на результат, так как распределение мощности по поверхности теплового потока имеет максимальные значенияв центре образца и убывает пропорционально радиусу. Разогрев образца происходит в центре пластинки и не достигает ее краев. Образец прогревается на глубину, не превышающую 1 мм. Таким образом основное внимание при расчете должно быть уделено граничным условиям на нагреваемой поверхности и на оси симметрии. Граничное условие на нагреваемой поверхности γ (второе из уравнений (2.2)) в цилиндрической геометрии принимает вид 1 ∂T 1 1 = ∂z 1z=0 W (t, r) - N (t, r) . λ(T ) Оно зависит от трех нелинейных коэффициентов. Напомним, что распределение мощности по поверхности теплового потока задается формулой (2.3), где A = 1/a2 = 0,03088523 мм-2 - постоянная, характеризующая величину радиуса пучка a, который имеет одно значение для каждой серии экспериментов и используется при расчете граничного условия для задачи (2.2). На каждом временном шаге численного моделирования значение переменной Wmax(t) берется из файла экспериментальных данных, индивидуального для каждого эксперимента. Потеря мощности на процесс испарения учитывается в виде 1 dm 1γ N (T 1 ) = Le S dt , где 1 dm - скорость испарения массы m на единицу площади поверхности S, а L - энтальпия S dt e парообразования. Потеря мощности рассчитывается по следующим формулам: Вт · мкс Le = 4,482 · 1012 кг , 1 dm ! M = P (T 1 ) = exp / 26,19104 - 83971,3 K \! 0,184 K кг · 10-12 . (3.2) S dt 1γ 1γ 2πRT 1 T 1γ 1γ 1 2π8,314T 1 мм2 · мкс 1γ Здесь T 1 1γ - температура на поверхности γ, R - универсальная газовая постоянная, а P (T 1 ) - давление насыщенного пара. Вывод уравнения (3.2) и обоснование пренебрежения излучением температуры поверхностью в энергетическом балансе более подробно представлены в работе [29]. Поскольку уравнение теплопроводности (3.1) традиционно аппроксимируется со вторым поряд- 1 ком точности по пространству, граничное условие на оси ∂T 1 ∂r 1r=0 = 0 в цилиндрической системе координат целесообразно аппроксимировать с тем же порядком точности следующим образом. Если значение во втором узле по радиусу разложить в ряд Тейлора Tn+1 n+1 1 ∂T 1n+1 h2 ∂ 2T 1 n+1 3 2k = T1k + h ∂r 1 + 1r=0 2 1 1 ∂r2 1r=0 + O(h ), то граничное условие на оси можно представить в виде: ∂T 1n+1 Tn+1 n+1 h ∂2T 1n+1 1 1 ∂r 1r=0 2k - T1k 1 1r=0 = h - 2 ∂r2 1 + O(h2). (3.3) Далее из уравнения (3.1) можно выразить вторую производную по радиусу: 1 ∂ r ∂r ( rλ(T ) ∂T ∂r ∂T = (c(T )ρ(T )+ Lmδ(T, ε)) ∂t ∂ - ∂z ( λ(T ) ∂T ∂z . (3.4) Расписав левую часть уравнения (3.4) в виде 1 ∂ r ∂r ( rλ(T ) ∂T = ∂r 1 ∂T λ(T ) 1. ∂r ∂ ( + λ(T ) ∂r ∂T ∂r , (3.5) 448 Д. Е. АПУШКИНСКАЯ, Г. Г. ЛАЗАРЕВА заметим, что вблизи оси симметрии решение практически постоянно, а теплопроводность λ(T ) ограничена. Поэтому применив к первому слагаемому в правой части (3.5) правило Лопиталя получим lim r→0 ( 1 λ(T ) r ∂T = ∂r ∂ ( λ(T ) ∂r ∂T , ∂r ∂ ( 2 λ(T ) ∂r ∂T ∂r ∂T = (c(T )ρ(T )+ Lmδ(T, ε)) ∂t ∂ - ∂z ( λ(T ) ∂T . ∂z Из последнего соотношения следует равенство ∂2T 1 ( ∂T ∂λ(T ) ∂T ∂ ( ∂T ∂r2 = 2λ(T ) (c(T )ρ(T )+ Lmδ(T, ε)) ∂t - 2 ∂r ∂r - ∂z λ(T ) , ∂z которое после подстановки в граничное условие (3.3) дает выражение 1 ( ∂λ(T ) ∂T 1n+1 Tn+1 - Tn+1 h 1 ∂T 1n+1 1 - 2 ∂r 1 = 2k 1k ∂r 1r=0 h - 4λ(T ) (c(T )ρ(T )+ Lmδ(T, ε)) 1 + ∂t 1r=0 ∂ ( + λ(T ) 1 ∂T 1n+1 + O(h2). ∂z ∂z 1 1r=0 Таким образом, итоговая система уравнений имеет следующий вид: T - T ⎧ n+1 n ⎪ ik ik 1 r Ln n+1 n+1 n n+1 n+1 = ⎪ τ r h2Cn ik i+1 k ik ik ik i i+1/2 + (T - T ) - ri-1/2L (T - T 1 k )l ⎪ i ik ⎪ 1 r ⎪ ( ) (T - - T )l , ⎪ + Ln Tn+1 n+1 - Ln n+1 n+1 C ⎪ ⎪ ⎪Cn n ik+1/2 ik n n i k+1 - Tik n ik-1/2 ik i k-1 ⎪ ik = c(Tik )ρ(Tik )+ Lmδ(Tik, ε), ⎪ ⎪Ln n ik = λ(Tik ), ⎪ ⎪ ⎪ Tn+1 4Tn+1 + 3T n+1 W n+1 n+1 i,3 - ⎪ i,2 i,1 i - Ni , = ⎪ 2h ⎪ ⎪ L n+1 i,1 n+1 - ⎨Wi = Wmax(t n+1 ) exp ( 2 Ari ), (3.6) n+1 ⎪ n+1 n+1 - i,N -1 T 4T + 3T ⎪ i,Nz -2 i,Nz z ⎪ = 0, ⎪ 2h ⎪ n+1 n+1 ( n+1 n ( ⎪ - T T ⎪ 2k 1k h Cn T2k - T2k 1 Ln n+1 n ⎪ h - 4Ln 1k 1k ⎪ ⎪ T - T - + ( ) τ h 1 k+1/2 i k+1 ik ⎪ n (Tn+1 n+1 )\l = 0, ⎪ ⎪ ⎪ n+1 ⎪ n+1 n+1 - L1 k-1/2 ik - Ti k-1 - T 4T + 3T ⎪ Nr -2,k Nr,k Nr -1,k ⎪ ⎪ ⎪⎩T 0 ik = T0. = 0, 2h Здесь индексы i и k принимают соответственно значения i = 2,... , Nr - 1 и k = 2,... , Nz - 1. Предложенная схема решения (3.6) имеет второй порядок аппроксимации по пространству и позволяет не выделять границу фазовых превращений. Расчет температуры во внутренних точках расчетной области реализуется с использованием схемы, стабилизирующей поправки, и метода прогонки [8], а коэффициенты рассчитываются по значениям температуры с предыдущего слоя. Уравнение теплопроводности (3.1) имеет дивергентный вид, что позволяет использовать консервативные схемы. Решение системы (3.6) с постоянными коэффициентами протестировано на известных аналитических решениях, с переменными коэффициентами протестировано на экспериментальных данных. При тестировании также учитывался известный факт, что разогрев образца вольфрама с плотностью мощности 103 Вт/мм2 в течение 1000 мкс приводит к разогреву образца до 2000 К. АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ СТЕФАНА 449 Предложенный алгоритм экономичен и прост в реализации, может быть использован для моделирования плавления других тугоплавких металлов. 4. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ РАСЧЕТОВ На основе схемы (3.6) производился расчет распределения температуры в поперечном сечении вольфрамовой пластинки. Таким образом учитывался прогрев вглубь материала, что и определяло температуру на поверхности наряду с падающей энергией и испарением, которые задаются граничным условием. Угловаясимметрия лазерного импульса определяла угловую симметрию распределения температуры. Мгновенный нагрев пластинки до 9000 K за счет воздействияна её поверхность теплового потока и последующее охлаждение сопровождались возникновением и исчезновением области расплава (см. рис. 3). РИС. 3. Изолинии температуры в поперечном сечении образца в процессе охлаждения в моменты времени 80 мкс (а) и 160 мкс (б), изолиния при температуре плавления выделена пунктиром, область расплава серым цветом. Область расплава в расчетах определяется по изолинии температуры, которая соответствует температуре плавления. На рис. 3 приведены изолинии температуры для типичного эксперимента. Изолиния, соответствующая температуре плавления, на рис. 3 выделена пунктиром, а область расплава отмечена серым цветом. Расчетная область на рисунке сильно растянута в осевом направлении, так как глубина прогрева мала по сравнению с радиусом расплава. На рис. 3а видно, что расплав в процессе нагрева в момент времени 80 мкс от момента начала воздействияимеет радиус около 4,5 мм и глубину около 0,02 мм. В процессе охлаждения(см. рис. 3б) область расплава сужается до 1,5 мм, но прогрев вглубь пластинки продолжается, глубина расплавленной области увеличивается в два раза к моменту времени 160 мкс. Сравнение других изолиний на рис. 3а и 3б показывает, что в два раза увеличиваетсяглубина разогретой области. Можно сделать вывод, что во время поверхностного охлаждения происходит разогрев вглубь пластинки, что означает смену знака скорости разогрева. Область расплава и температуру поверхности возможно фиксировать в ходе экспериментов. Например, анализ влияния испарения на процесс нагрева был проведен на основе измерений радиуса расплавленной области [10]. Результаты расчетов показали полное совпадение динамики расчетного и экспериментального радиуса расплава в случае учета испаренияна границе пластинки вольфрама. Из экспериментальных данных известно, что за счет постоянного процесса испарения образец не разогреваетсявыше 9000 K. Процесс испарениявносит свой вклад в разогрев вещества начиная с 2000 K, то есть при температурах много ниже температуры плавления. На рис. 4а представлена временная зависимость максимума плотности мощности нагрева электронным пучком, которая является входными данными для расчета выстрелов 679, 680 и 681. На рис. 4б представлены графики температуры поверхности для этих экспериментов. Моменты, для которых выводятся результаты расчета, соответствуют времени измерения. Эти моменты времени на графике 4а обозначены заштрихованными прямоугольниками, ширина основания которых равна времени экспозиции измерительного прибора. Графики используются при анализе измерений и планировании экспериментов. 450 Д. Е. АПУШКИНСКАЯ, Г. Г. ЛАЗАРЕВА РИС. 4. Динамика плотности мощности в центре зоны воздействияи моменты измерения(а), графики температуры поверхности (б). Данные приведены длявыстрелов 679, 680 и 681. Результаты расчета радиального распределения температуры по поверхности образца (рис. 4) полностью согласуются с экспериментальными данными и аналитическими оценками из работы [10]. Разработка дискретной модели позволяет проводить расчеты для экспериментов с большим временем задержки между измерениями и большим временем экспозиции. На рис. 4а показаны моменты измерения 100-110 мкс, 500-550 мкс, 1500-1800 мкс, 5000-6000 мкс для со временем экспозиции 10, 50, 300 и 1000 мкс. То есть первый калибровочный кадр был сделан в интервале от 100 до 110 мкс со временем экспозиции 10 мкс. Графики температуры поверхности были рассчитаны для моментов времени 105 мкс, 525 мкс, 1650 мкс и 5500 мкс (рис. 4б). Например, первый калибровочный кадр был сделан в интервале от 100 до 110 мкс со временем экспозиции 10 мкс, соответствующий ему расчет был сделан для момента времени 105 мкс. Продолжительность нагрева составила 124 мкс, 132 мкс и 130 мкс при воздействиях 679, 680 и 681 соответственно. Небольшие отличия результатов обусловлены разной энергией и длительностью воздействия. В результате нагрева температура образцов приближается к 8000 К. После окончаниявоздействия материал остывает, температура поверхности снижаетсядо 4000 K, глубина распространения тепла в объеме образца составляет менее четверти миллиметра. На графике прямая линия в точке плавления показывает площадь расплава. Из многочисленных контрпримеров хорошо известно, что вторые производные от температуры T по пространственным переменным r и z, так же как и первая производная от T по времени, могут иметь скачок на свободной границе. Подчеркнем, что графики на рис. 4б как раз показывают изменение характера температурной кривой в точках раздела расплава и твердой фазы, что является следствием свойств решения задачи Стефана (2.2). Спад ожидаемой кривизны обусловлен в первую очередь потерями энергии на испарение. 5. ВЫВОДЫ Представленное численное решение задачи Стефана для расчета температуры образца вольфрама, нагреваемого импульсным электронным пучком, показывает хорошее соответствие как с результатами натурного эксперимента, так и с аналитическими фактами теории задач со свободными границами. Разработанный алгоритм решениязадачи Стефана в дальнейшем послужит основой для расчета тока в образце, который рассматривается как возможный источник наблюдаемого в эксперименте вращения вещества. Предложенный алгоритм отличается простотой и экономичностью. Он может быть применен к моделированию процесса плавления других тугоплавких металлов.

Об авторах

Д. Е. Апушкинская

Российский университет дружбы народов

Автор, ответственный за переписку.
Email: apushkinskaya_de@pfur.ru
Москва, Россия

Г. Г. Лазарева

Российский университет дружбы народов

Email: lazareva_gg@pfur.ru
Москва, Россия

Список литературы