自由边界问题

自由边界问题（free-boundary problem）是一类需要同时求出区域内部的未知函数和区域边界位置的问题。边界通常由解本身决定，而非预先给定。冰块融化时的固液界面、液体的自由表面和某些接触问题中的接触区域，都可以产生这样的数学结构。

以从一侧加热的固体为例：即使知道热量怎样在液体中传导，也还不知道已经融化了多厚；而融化区域的厚度又决定温度方程在哪个区间求解。因此温度和界面必须联立计算。这个例子称为 Stefan 问题，下面推导一个能完整求出的理想化情形。

未知边界与给定运动边界

设液相占据 $0<x<s(t)$。若 $s(t)$ 已由实验装置或几何运动指定，温度方程是在给定运动区域上的问题；若 $s(t)$ 也由热平衡决定，它就是自由边界。这里“自由”指位置待求，并不表示界面没有约束。

一个完整模型通常包含区域内的微分方程、固定边界上的条件、自由界面上的条件和初始状态。只给出温度方程及两端温度，不能额外确定任意运动的 $s(t)$；还需要一条控制界面运动的规律。相变中的这条规律是潜热收支。Stefan 问题的单相与双相设置及相似解参数可参见 RWTH Aachen 的课程资料。

一相融化模型

考虑半无限固体 $x>0$，初始温度恰好为熔点 $T_m$。从 $t>0$ 起，左端壁面保持温度 $T_h>T_m$，液相在左侧形成，界面向右推进。假设材料为纯物质、界面尖锐，忽略液体对流及固液密度差，物性参数为常数。固体始终保持 $T_m$，只需求液体中的温度。

记密度为 $\rho$，比热容为 $c$，导热系数为 $\kappa$，单位质量熔化潜热为 $\ell >0$，热扩散率为 $\alpha =\kappa /(\rho c)$。液相温度满足热方程 $$T_t=\alpha T_{xx},0<x<s(t),t>0,$$ 边界条件为 $$T(0,t)=T_h,T(s(t),t)=T_m,s(0)=0.$$ 导热通量是 $q=-\kappa T_x$。在时间 $\mathrm{d}t$ 内，单位横截面积新融化的质量为 $\rho s^{\prime }(t)\mathrm{d}t$，需要潜热 $\rho \ell s^{\prime }(t)\mathrm{d}t$。固相没有温度梯度，故界面收到的液相热通量全部用于融化： $$\rho \ell s^{\prime }(t)=-\kappa T_x(s(t{)}^{-},t).$$ 上标负号表示从液相一侧取极限。温度沿正方向降低，$T_x<0$，于是右端为正，符合界面向右运动。若两相都有温度梯度，就要使用两侧通量之差，而不能只保留这一侧。

相似变量与温度剖面的推导

扩散的自然长度尺度为 $\sqrt{\alpha t}$。因此尝试 $$\eta =\frac{x}{2\sqrt{\alpha t}},s(t)=2\lambda \sqrt{\alpha t},T(x,t)=T_m+\Delta TF(\eta ),\Delta T=T_h-T_m.$$ 这里 $\lambda >0$ 仍是待求参数。计算导数，得到 $$T_t=-\frac{\Delta T}{2t}\eta F^{\prime }(\eta ),\alpha T_{xx}=\frac{\Delta T}{4t}{F}^{″}(\eta ).$$ 所以热方程化为 $F^{″}+2\eta F^{\prime }=0$。令 $g=F^{\prime }$，则 $g^{\prime }+2\eta g=0$。乘以积分因子 $e^{{\eta }^2}$ 后，$(e^{{\eta }^2}g{)}^{\prime }=0$，故 $g=C{e}^{-{\eta }^2}$。定义误差函数 $$\mathrm{erf}(z)=\frac{2}{\sqrt{\pi }}{\int }_0^z{e}^{-r^2}\mathrm{d}r.$$ 再用 $F(0)=1$、$F(\lambda )=0$，便得到 $$F(\eta )=1-\frac{\mathrm{erf}(\eta )}{\mathrm{erf}(\lambda )},$$ $$T(x,t)=T_h-\Delta T\frac{\mathrm{erf}(x/(2\sqrt{\alpha t}))}{\mathrm{erf}(\lambda )},0\le x\le s(t),t>0.$$ 分母保证界面温度恰好等于 $T_m$。该表达式只用于液相；界面右侧的固相温度仍为 $T_m$。

纵轴为归一化超温 θ=(T−Tm)/(Th−Tm)，取无量纲 α=1、Ste=1。t=0.25、1、4 的剖面在各自界面处到达熔点；界面右侧是等温固体。

由潜热确定界面速度

相似形式还给出 $$s^{\prime }(t)=\frac{\lambda \sqrt{\alpha }}{\sqrt{t}},T_x(s(t{)}^{-},t)=-\frac{\Delta T{e}^{-{\lambda }^2}}{\mathrm{erf}(\lambda )\sqrt{\pi \alpha t}}.$$ 将两式代入界面热平衡，并使用 $\kappa =\rho c\alpha$，消去 $\rho \sqrt{\alpha }/\sqrt{t}$ 后得到 $$\lambda e^{{\lambda }^2}\mathrm{erf}(\lambda )=\frac{\mathrm{S}\mathrm{t}\mathrm{e}}{\sqrt{\pi }},\mathrm{S}\mathrm{t}\mathrm{e}=\frac{c\Delta T}{\ell }.$$ $\mathrm{S}\mathrm{t}\mathrm{e}$ 是 Stefan 数，比较升温所对应的显热尺度与熔化潜热。它没有量纲。左侧函数记为 $G(\lambda )$，求导得 $$G^{\prime }(\lambda )=e^{{\lambda }^2}(1+2{\lambda }^2)\mathrm{erf}(\lambda )+\frac{2\lambda }{\sqrt{\pi }}>0(\lambda >0).$$ 又有 $G(0)=0$、$G(\lambda )\to \mathrm{\infty }$，所以每个正 Stefan 数恰有一个正相似参数。这证明了本相似解族中的参数唯一性，并非对所有自由边界模型作一般的唯一性断言。

在教学用无量纲参数 $\alpha =1,\mathrm{S}\mathrm{t}\mathrm{e}=1$ 下，对单调方程用二分法求根，得到 $$\lambda \approx 0.6200626333,s(t)\approx 1.240125267\sqrt{t}.$$ 因此 $t=0.25,1,4$ 时的界面位置约为 $0.6201,1.2401,2.4803$。时间增大四倍，融化厚度增大两倍；界面速度反而降低，因为热量需要穿过更厚的液层。

s(t)=2λ√t，标出三个与温度剖面相同的时刻。曲线在初始时刻的无界斜率来自瞬间升高壁温、零初始液层的理想化。

一个可检验的变式是：保持物性及壁温不变，潜热加倍。此时 Stefan 数减半，因 $G$ 严格递增，新的 $\lambda$ 必较小，界面在相同时间推进得更少。这一步无需先求出新根；它直接来自模型中的热量收支与参数方程。

总热量的交叉检验

相似解还应满足整个融化区的热量收支。相对初始固体，单位面积储存的能量为 $$E(t)=\rho \ell s(t)+\rho c{\int }_0^{s(t)}(T-T_m)\mathrm{d}x.$$ 第一项是潜热，第二项是液体高于熔点的显热。对移动上限求导，因界面处 $T-T_m=0$，积分的上限项为零；再代入热方程，得到 $$E^{\prime }=\rho \ell s^{\prime }+\kappa {\int }_0^{s(t)}{T}_{xx}\mathrm{d}x=\rho \ell s^{\prime }+\kappa T_x(s^{-})-\kappa T_x(0^{+})=-\kappa T_x(0^{+}).$$ 界面处的两项由 Stefan 条件抵消，只剩壁面输入热通量。这一计算能检查符号，也能检查数值程序是否把潜热和显热重复计入。实际数值模型常以尖锐界面相似解作为相变算法的测试基准；将相变摊在一段温度区间内时，还需考察该区间及网格、时间步的影响。

固定参考区间上的计算

为在固定网格上计算，可令 $\xi =x/s(t)$、$v(\xi ,t)=T(s(t)\xi ,t)$。液相区间变为 $0<\xi <1$。固定物理位置求时间导数时，${\xi }_t=-\xi s^{\prime }/s$，因此 $$T_t=v_t-\frac{\xi s^{\prime }}{s}{v}_{\xi },T_{xx}=\frac{1}{s^2}{v}_{\xi \xi }.$$ 变换后的方程及界面条件为 $$v_t=\frac{\alpha }{s^2}{v}_{\xi \xi }+\frac{\xi s^{\prime }}{s}{v}_{\xi },\rho \ell s^{\prime }=-\frac{\kappa }{s}{v}_{\xi }(1,t).$$ 新增的一阶项来自坐标移动，不是原材料中另加了一股流。$s(t)$ 仍须求解；变换只是固定计算区域，并没有消除自由边界的未知性。

参考点 ξ=0、1/4、1/2、3/4、1 保持不变，对应物理点 x=s(t)ξ 随界面伸展。图中两行分别为 t=1 与 t=4。

由于本例 $s(0)=0$，该坐标变换在零时刻没有定义。数值计算可从已知相似解的某个正时刻初始化，或专门处理起始层，不能把含 $1/s$ 的式子直接代入零时刻。若存在对流、材料物性变化、非平衡界面温度或复杂几何，平方根规律也未必继续成立。