热源周期振荡条件下一维融化问题的数值研究

曲良辉，杨金根，邢 琳，令 锋，董 丽

(1. 中原工学院 理学院，河南 郑州 450007；2. 信阳师范学院 数学与信息科学学院，河南 信阳 464000；3. 肇庆学院 数学与统计学院，广东 肇庆 526061)

0 引言

在自然界和各种工业生产中，经常会遇到这些过程，如冰的融化、金属的重结晶、合金的凝固、食物的冷藏、液滴的蒸发和粒子的扩散等，它们都涉及相变传热问题[1-3]．其特点是区域内存在一个随时间变化的两相界面，在该界面上吸收或放出热量．这类问题通常被称为Stefan问题，也称为移动边界问题，因为两相界面的位置是移动的．在数学上相变传热问题是一个强非线性问题，由于两相界面的位置未知，因此界面上能量守恒条件也是非线性的．求解这类问题时，一般多用近似方法或数值方法求解[4-6]．数值方法求解时，通常应用空间坐标变换，把移动区域模型转化为固定区域模型进行间接求解[7-8]．

本文以一维介质的融化为模型，考虑在初始位置x=0处有一周期振荡的热源，融化从x=0处开始，固液交界面的移动边界X(t)沿x正方向移动．令固态区域总保持恒温并且恰好为相变温度．而液态区域，通过空间坐标变换，将移动边界转化为固定边界，采用显式和隐式两种有限差分格式数值模拟在融化过程中移动边界的运动和液态介质内温度场的分布，进而比较两种有限差分方法的数值稳定性、计算量和计算效率，为该类型的相变传热问题提供参考．

1 数学模型

令x=x′/X′(t′)，一维移动边界问题能够通过空间坐标的简单拉伸转化为固定边界问题[7]．这里x′是空间变量,0≤x′≤X′(t′);X′(t′)是移动边界位置;x是无量纲化的空间变量,0≤x≤1．

引进无量纲变量和参数，在x=0处温度随时间周期振荡的一维融化问题，其固定区域无量纲化模型可表示为[8]：

(1)

(2)

边界和初始条件为

T(x=1,t)=0,t>0,

(3)

T(x=0,t)=1+εsin(ωt),t≥0,

(4)

X=0,t≤0,

(5)

2 两种有限差分格式的建立

对时间进行离散时，先取定一个小的时间t0作为计算的初始时刻，令tm=t0+mΔt，其中Δt为时间步长，m=0,1,2,…．对空间区域x∈[0,1]进行N等分，令xj=jh，j=0,1,2,…,N，其中h为空间步长，并且x0=0，xN=1．

2.1 显式有限差分格式的建立

微分方程组(1)-(5)的显式有限差分格式为[8]：

j=1,2,…,N-1,

(6)

m=0,1,2,…,

(7)

(8)

(9)

X0=X(t0).

(10)

其截断误差为O(Δt+h2)．应用冻结系数和Fourier稳定性判定方法，上述显格式在满足如下条件时是数值稳定的[9]：

2.2 隐式有限差分格式的建立

微分方程组(1)-(5)的隐式有限差分格式为：

j=1,2,…,N-1,

(11)

m=0,1,2,…,

(12)

(13)

(14)

X0=X(t0).

(15)

其截断误差为O(Δt+h2)．借助冻结系数的方法，可以证明该隐格式具有绝对的数值稳定性[9]．

3 数值结果和讨论

这里采用显式和隐式两种有限差分格式,对x=0处温度周期振荡(振荡振幅ε=0.5，振荡频率ω=π/2)的条件下介质融化过程中移动边界的运动和液态介质内温度场的分布进行数值模拟计算，并对两种有限差分格式的数值稳定性、计算量和计算效率进行比较．这里所进行的数值实验是用MATLAB 7.1编程、在系统为Microsoft Windows XP，CPU为Pentium(R)4/2.80 GHz，内存为1 G的个人计算机上进行的．

考虑纯物质的融化问题，对于热源为定常温度T(x=0,t)=1 (t>0)(即ε=0)的这种情况，微分方程组(1)-(5)的精确解为[10]：

(16)

(17)

其中λ满足下面的超越方程

(18)

对于热源周期振荡的一维融化问题，需要说明的是为了避免t=0时的特殊情况，取t=0之后的某个时刻t0作为计算过程的初始时刻，应用式(16)和(17)近似地对时刻t0处的温度场和移动边界位置进行初始化，这样显格式和隐格式两种有限差分格式的数值计算过程就能顺利向前推进．

由显式有限差分格式稳定性的条件可知，显格式对时间步长的选取是有限制的且要求时间步长非常小．这里取初始时刻为t0=0.01，空间步长为h=1/N=0.1，振荡振幅为ε=0.5.通过数值实验可知，当时间步长取Δt=0.000 02时，对于Stefan数分别取Ste=0.2、1.0和2.0，都可以保证显式差分格式的数值稳定性．若时间步长取Δt=0.000 1时，对于Ste=2.0这种情况是符合上面所得显格式稳定性的充分条件，而对于Ste=1.0这种情况尽管不符合上面所给显格式稳定性的充分条件，但通过数值计算发现还是可以保证显格式的数值稳定性．这说明了上述所给显格式的稳定性条件是充分的，没有包括时间步长Δt的所有可能取值，主要是由于稳定性分析时采用了“冻结系数”进行预处理的结果．而隐式有限差分格式是一种绝对稳定的差分格式[9]，对时间步长和空间步长的选取是没有限制的．图1显示了在Ste=1.0、时间步长Δt分别取0.1、0.01和0.001时应用隐式有限差分格式所获得的移动边界位置随时间变化的关系曲线．可以看出所得的三条曲线近似程度非常高.特别地，当时间步长为0.01和0.001时，所对应的两条曲线完全重合；另一方面从数值上也进一步表明了所构造的隐式有限差分格式的绝对稳定性．

图2数值模拟了显式差分格式和隐式差分格式对三个不同的Stefan数(Ste=0.2、1.0和2.0)所获得的融化介质移动边界的位置X(t)随时间t的变化关系．从图2(a)和图2(b)可以看出，随着Stefan数的增加，移动边界的运动速度也随着增大，而且在x=0处周期振荡的热源对介质相变速度的影响会在相变初期较短的时间段内比较明显．随着时间的增加，在一段时间之后，移动边界的位置变化逐渐与时间的平方根成正比．此外，从图2中还可以发现，使用这两种差分格式所得的数值结果吻合地非常好，都能够很好地反映在热源周期振荡条件下移动边界的位置随时间的变化规律．

图1 隐格式下不同时间步长所对应的移动边界的运动

图2 不同Stefan数下两种差分方法得到的移动边界的位置与时间的关系

图3数值模拟了应用两种差分格式所得到的在不同时刻液态介质内的温度分布．从图3中可以看出，在x=0处周期振荡的温度强烈影响着液态介质内的温度分布．通过比较图3(a)和图3(b)可以看出，这两种差分格式所模拟的介质内温度场的分布高度吻合．结合图2和图3，可推断出采用这两种有限差分格式对热源周期振荡条件下融化问题进行数值计算都是可行的．

图3 两种差分方法得到的液态介质内不同时刻的温度分布

下面从计算量和计算效率上对这两种有限差分格式进行比较．表1显示的是h=0.1和Ste=1.0时两种有限差分格式在计算到t=25.0时所需的计算量和计算时间．可以看出，对于显格式，如果采用时间步长为Δt=0.000 1，则计算到t=25.0这一时刻需要循环计算249 900次，计算量非常大；如果取时间步长为Δt=0.001时进行计算，数值结果不再收敛，出现了不稳定的情况．然而对于绝对稳定的隐格式来说，可以采用比较大的时间步长．如果采用时间步长为Δt=0.01，计算到t=25.0这一时刻只需要循环计算2 499次；若采用时间步长为Δt=0.1进行计算，仍然可以保证数值结果的稳定性和收敛性，这时只需要循环计算250次，大大减少了计算量．此外，从表1中还可以看出，仅仅计算到t=25.0这一时刻，显格式所需的运行时间竟然达到13 303.137 805 s.而隐式差分格式在保证稳定性和收敛性的前提下，时间步长可以取Δt=0.1进行计算，虽然每一次的循环计算都需要反复迭代才能获得满足误差要求的结果，但只需5.242 976 s即可获得t=25.0这一时刻的结果．因此，相比于显式差分格式，隐式差分格式在介质融化的理论分析中可以选取比较大的时间步长进行计算，能够大大减少计算量和计算时间而更具有应用价值．

表1 h=0.1和Ste=1.0时两种差分格式计算量和计算时间的对比

4 结论

通过空间坐标变换，把移动区域模型转化为固定区域模型，分别采用显式有限差分格式和隐式有限差分格式求解了热源周期振荡条件下的一维融化问题．通过对融化过程中移动边界的运动及液态介质内温度场分布的数值模拟发现这两种差分格式的数值结果高度吻合，但是在数值稳定性、计算量和计算效率方面隐式差分格式明显优于显式差分格式．