曲良辉 邢 琳 钱德亮 令 锋

低温边界条件下柱型凝固问题的数值研究

曲良辉 邢 琳 钱德亮 令 锋

通过考虑移动边界每次向外推进固定距离时所需要的时间，构造定空间步长方法求解了两种不同低温边界条件下的向外柱型凝固问题，数值模拟了相变过程中移动边界的运动及介质内温度场的变化。数值结果表明采用定空间步长方法求解柱型凝固问题是数值稳定的，并且具有较高的精度，对于Stefan问题的研究具有重要的参考价值。

Stefan问题也称为移动边界问题，其特点是区域内存在着一个随时间运动的边界，如冰的融化、金属铸件的凝固、微粒扩散和种群繁衍等。求解Stefan问题，由于问题的非线性性质和移动边界的未知位置，只有一些特殊情况才能获得解析解，通常只能用近似方法或数值方法求解。对于柱状模型下的相变传热问题，1981年Lunardini借助热平衡积分技巧，应用等效热扩散率方法获得了具有定常热源边界条件的柱型相变问题的一个近似解析解。2009年，Caldwell等应用焓法、边界固定法、摄动法和热平衡积分法分别求解了柱型相变传热问题，并就定常边界条件的情况，对数值结果进行了分析和比较。

本文考虑在低温边界条件下，通过计算移动边界每次向外推进固定距离时所需要的时间，构造了一种定空间步长方法数值求解向外柱型凝固问题，并对移动边界的运动和固态介质内的温度场进行了数值模拟，分析了低温边界条件对相变区域内的温度变化及移动边界运动的影响。

数学模型

考虑初始温度恰好为相变温度Tf的某液态介质由于边界处温度的变化所导致的向外柱型凝固问题（图1），在相变过程中体积与物性均不发生变化的情况下，该问题的数学模型为

满足的条件为

引进无量纲量

定空间步长方法

记Δr 为固定的空间网格步长，tj（t0=0）是对应于移动边界位于Rj=1+jΔ r （j=0，1，…）处的时间。考虑相变界面移动到RN=1+NΔ r 处所对应的时间tN，此时各节点rj=1+jΔ r 处的固相介质的温度记为由于控制方程可以简化为

因此，应用向后差分格式离散时间导数项和中心差分格式离散空间导数项，则控制方程的离散形式为

Stefan条件（7）的离散形式为

引入两个变量a=1Δr 和b=ΔrΔtN-1，式（11）和（12）分别变为

边界条件（8）-（10）的离散形式分别为

数值结果与讨论

图1　向外柱型凝固示意图

在Stefan数Ste=0.2和f（ t）=0的情况下，应用定空间步长方法，图2（a）数值模拟了移动边界位置随时间变化的运动曲线和图2（b）数值模拟了当移动边界分别向外移动到R（ t）=1.2，1.5和1.8时，相变区域内温度的变化曲线，其中空间变量步长固定为Δr=0.001。表1对定空间步长方法所得到的不同时刻处移动边界的位置与文献中的相应结果进行了比较。从表中可以看出，在Stefan数Ste=0.2和边界条件f（ t）=0的情况下，定空间步长方法所得的结果与文献中的焓法和摄动法所得到的相应结果非常接近，其绝对误差均不超过0.001，从而表明了该数值方法求解柱型凝固问题的可行性和有效性。

图2　在Ste=0.2和f（ t）=0的情况下移动边界的运动和相变区域内的温度分布

表1　在Ste=0.2和f（ t）=0的情况下移动边界位置的比较

在Stefan数Ste=0.2和f（ t）=-t2的情况下，图3（a）数值模拟了应用定空间步长方法所得到的移动边界位置随时间变化的运动曲线，图3（b）数值模拟了当移动边界分别向外移动到R（ t ）1.2，1.6和2.0时，所得到的相变区域内温度的变化曲线，其中使用的空间变量步长为Δr=0.001。从图3（a）中可以看出，在一段时间之后，移动边界的位置随时间的变化几乎保持线性的运动趋势，也就是说移动边界的运动速度几乎不变。从图3（b）中可以看出，相变区域内的温度场在快速变化，当移动边界的位置从R（ t）=1.6向外推进到R（ t）=2.0时，相变区域内的温差从大约1.8增加到了约4.3，这是受低温边界条件影响的结果。

图3　在Ste=0.2和f（ t）=-t2的情况下移动边界的运动和相变区域内的温度分布

结束语

通过考虑相变界面每次向外推进固定距离时所需要的时间，构造了定空间步长方法用来求解向外柱型凝固问题，并模拟了相变过程中移动边界的运动和固态介质内温度场的变化。通过与已有结果比较可知，应用定空间步长方法求解向外凝固问题是可行的，并且具有较高的精度，从而表明了定空间步长方法在求解Stefan问题时所具有的重要应用价值。

DOI：10.3969/j.issn.1001-8972.2015.15.005