Задача Стефана

Рассматривается полубесконечный одномерный образец льда, находящийся при температуре плавления ${\displaystyle u(x,t)\equiv 0}$ для ${\displaystyle x\in [0,+\infty )}$. Положение границы между льдом и водой обозначается через ${\displaystyle s(t)}$. Поток тепла ${\displaystyle f(t)}$ действует на левой границе, что приводит к плавлению льда и увеличению области ${\displaystyle [0,s(t)]}$, занимаемой водой. В безразмерных переменных постановка имеет вид^[3][4]:

• ${\displaystyle {\frac {\partial u}{\partial t}}={\frac {\partial ^{2}u}{\partial x^{2}}},\quad 0<x<s(t),\;t>0}$ — уравнение теплопроводности, описывающее изменение температуры;
• ${\displaystyle -{\frac {\partial u}{\partial x}}(0,t)=f(t),\;t>0}$ — условие Неймана на левом конце области, задающее поток тепла на входе;
• ${\displaystyle u{\big (}s(t),t{\big )}=0,\;t>0}$ — условие Дирихле на межфазной границе вода — лёд;
• ${\displaystyle {\frac {\mathrm {d} s}{\mathrm {d} t}}=-{\frac {\partial u}{\partial x}}{\big (}s(t),t{\big )},\;t>0}$ — условие Стефана, определяющее скорость движения межфазной границы;
• ${\displaystyle u(x,0)=0,\;x\geq 0,\qquad s(0)=0}$ — начальные условия.

Задача называется однофазной, поскольку температура в твёрдой фазе тождественно равна температуре плавления и уравнение теплопроводности решается только в жидкой области.

В размерных переменных для двухфазной задачи о плавлении в каждой из фаз ${\displaystyle i=1,2}$ решается уравнение:

${\displaystyle \rho _{i}c_{i}{\frac {\partial T_{i}}{\partial t}}={\frac {\partial }{\partial x}}\left(\lambda _{i}{\frac {\partial T_{i}}{\partial x}}\right),}$

а на границе ${\displaystyle x=s(t)}$ выполняются условия ${\displaystyle T_{1}=T_{2}=T_{m}}$ и условие теплового баланса:

${\displaystyle \rho L{\frac {\mathrm {d} s}{\mathrm {d} t}}=\lambda _{1}{\frac {\partial T_{1}}{\partial x}}{\Big |}_{s(t)-0}-\lambda _{2}{\frac {\partial T_{2}}{\partial x}}{\Big |}_{s(t)+0},}$

где ${\displaystyle \rho }$ — плотность, ${\displaystyle c_{i}}$ — удельная теплоёмкость, ${\displaystyle \lambda _{i}}$ — коэффициент теплопроводности, ${\displaystyle L}$ — удельная теплота фазового перехода, ${\displaystyle T_{m}}$ — температура плавления.

Для анализа задачи вводятся ключевые безразмерные параметры^[6]:

• число Стефана (${\displaystyle \mathrm {Ste} ={\frac {c\,\Delta T}{L}}}$): характеризует отношение явной (чувствительной) теплоты к скрытой теплоте фазового перехода. При ${\displaystyle \mathrm {Ste} \ll 1}$ применима квазистационарная аппроксимация;
• число Фурье (${\displaystyle \mathrm {Fo} ={\frac {at}{L_{c}^{2}}}}$): безразмерное время, характеризующее отношение скорости теплопроводности к скорости накопления тепла;
• число Био (${\displaystyle \mathrm {Bi} ={\frac {hL_{c}}{\lambda }}}$): определяет соотношение между внутренним термическим сопротивлением тела и сопротивлением теплоотдаче на его поверхности.

Рассматривается процесс диффузионного взаимодействия в бинарной металлической системе ${\displaystyle A-B}$ с ${\displaystyle \alpha }$- и ${\displaystyle \beta }$-фазами, представляющими собой регулярные твёрдые растворы. Через ${\displaystyle s(t)}$ обозначено положение подвижной межфазной границы, тогда ${\displaystyle \alpha }$-фаза занимает область ${\displaystyle 0\leq x<s(t)}$, а ${\displaystyle \beta }$-фаза — область ${\displaystyle s(t)<x\leq l}$^[21].

• ${\displaystyle {\frac {\partial N_{B}}{\partial t}}=D^{\alpha }\cdot {\frac {\partial ^{2}N_{B}}{\partial x^{2}}}}$ — уравнение, описывающее изменение концентрации в ${\displaystyle \alpha }$-фазе;
• ${\displaystyle {\frac {\partial N_{B}}{\partial t}}=D^{\beta }\cdot {\frac {\partial ^{2}N_{B}}{\partial x^{2}}}}$ — уравнение, описывающее изменение концентрации в ${\displaystyle \beta }$-фазе;
• ${\displaystyle \left(N_{B}^{\beta }-N_{B}^{\alpha }\right)\cdot {\frac {\mathrm {d} s}{\mathrm {d} t}}=-D^{\beta }\cdot {\frac {\partial N_{B}}{\partial x}}{\Bigr |}_{s(t)+0}+D^{\alpha }\cdot {\frac {\partial N_{B}}{\partial x}}{\Bigr |}_{s(t)-0}}$ — условие Стефана, определяющее скорость движения межфазной границы;
• ${\displaystyle {\frac {\partial N_{B}}{\partial x}}{\Bigr |}_{x=0}=0,\quad {\frac {\partial N_{B}}{\partial x}}{\Bigr |}_{x=l}=0}$ — граничные условия отсутствия потока вещества через внешние границы,

где ${\displaystyle N_{B}(x,t)}$ — концентрация атомов сорта ${\displaystyle B}$, ${\displaystyle D^{\alpha }}$ и ${\displaystyle D^{\beta }}$ — коэффициенты диффузии в фазах. Значения ${\displaystyle N_{B}^{\alpha }}$ и ${\displaystyle N_{B}^{\beta }}$ определяются диаграммой состояния системы и в предположении локального равновесия на границе считаются постоянными.

Исходная задача записывается в обобщённой формулировке в виде единого уравнения с разрывными коэффициентами (например, энтальпийная формулировка). Для построения численного алгоритма проводят процедуру сглаживания разрывных коэффициентов. Подход был предложен А. А. Самарским и Б. М. Будаком^[19][20]. Недостатком является «размытие» фронта на несколько ячеек сетки.

В этих методах положение границы отслеживается явно. К ним относятся:

• метод выпрямления фронта (front-fixing method): заменой переменных ${\displaystyle \xi =x/s(t)}$ подвижная область отображается в фиксированную, но уравнения становятся нелинейными^[17];
• метод ловли фронта в узел сетки: шаг по времени подбирается так, чтобы за него фронт сместился ровно на один пространственный шаг^[18];
• метод динамически адаптирующихся сеток: узлы расчётной сетки перемещаются вместе с фронтом, обеспечивая высокое разрешение в области фазового перехода^[24].

Современный альтернативный подход, в котором резкая граница раздела фаз заменяется тонким переходным слоем с непрерывным параметром порядка ${\displaystyle \phi }$ (где ${\displaystyle \phi =1}$ в одной фазе, ${\displaystyle \phi =-1}$ в другой). Эволюция ${\displaystyle \phi }$ описывается уравнениями типа Аллена — Кана или моделью Кобаяши, сопряжёнными с уравнением теплопроводности. Этот метод естественным образом описывает сложные топологические изменения (слияние, разрыв границ) и дендритный рост без необходимости явного отслеживания координат фронта^[23].

Для водоёма при температуре замерзания и температуре воздуха ${\displaystyle T_{a}=-10\,^{\circ }\mathrm {C} }$, при ${\displaystyle \lambda =2{,}22\ {\text{Вт/(м·К)}}}$, ${\displaystyle \rho =917\ {\text{кг/м}}^{3}}$ и ${\displaystyle L=3{,}34\cdot 10^{5}\ {\text{Дж/кг}}}$, формула Стефана даёт толщину льда около 0,21 м за сутки и около 0,56 м за неделю. Поскольку упрощённая модель пренебрегает теплоотдачей от воды и теплоёмкостью уже образовавшегося льда, она даёт оценку толщины льда сверху^[6].

При контакте двух металлов на границе образуется слой интерметаллида. В предположении диффузионного контроля роста и постоянных равновесных концентраций на границах толщина слоя изменяется по параболическому закону ${\displaystyle \Delta s\propto {\sqrt {t}}}$, а константа параболического роста выражается через коэффициенты диффузии и скачки концентраций, определяемые диаграммой состояния.