方腔内石蜡固液相变及其非线性特性分析

张 政，杨 茉，李易蓉

(上海理工大学 能源与动力工程学院，上海 200093)

相变储能技术通过余热回收利用来提高能源的利用率，进而可解决能源的浪费问题[1-3]。针对相变传热，国内外学者进行了很多研究。杜雁霞等[4]针对方腔内相变材料融化过程中自然对流的影响进行了实验研究。邹得球等[5]对石蜡相变过程进行了数值模拟，研究了相变的4个阶段对温度的影响。郭茶秀等[6]对充满方腔的相变材料进行了数值模拟，发现在两侧面和底面进行加热时融化速度最快。Zennouhi等[7]研究了在不同倾角下方腔内相变材料金属镓的传热过程和融化现象。Kousksou等[8]研究了温度均匀的垂直波浪形表面的熔化过程。战乃岩等[9]通过数值模拟方法讨论了方腔内流体流动、换热的静态分叉和振荡等非线性现象。Chen等[10]将方腔内的填充石蜡作为材料进行数值模拟，采用界面追踪法对以对流控制为主导的融化进行分析。Shamsundar等[11]采用焓模型进行了数值模拟，结果表明焓模型的数学表达式等价于固液两相区和固液界面的常规守恒方程。Korti等[12]研究发现倾斜角对相变材料融化过程中的自然对流有较大影响。

目前，关于相变材料融化流场中非线性特性的研究很少，而不同的融化结果会影响达到融化稳定阶段所需的时间。笔者将蓄热器抽象成一个底部加热、顶部冷却的方腔，采用焓法研究了融化过程中自然对流对融化的影响，通过保持上下壁面温度恒定，仅改变模型高度来改变瑞利数Ra，研究流场中速度的变化规律，并证明在相变过程中其数值结果存在多解等非线性特性。

1 模型的建立

1.1 模型

如图1所示，设置方腔模型长B为20 mm，高度A分别为5 mm、9.2 mm和20 mm。下壁面为加热端，温度为373.15 K；上壁面为冷却端，温度为335.15 K；左右壁面绝热，壁面厚度忽略不计。方腔内充满的相变材料为石蜡，模型整体水平放置。采用Fluent软件进行数值模拟，模型为二维非稳态，采用基于焓法的Solidification/Melting模型进行数值计算。采用SIMPLEC算法求解速度与压力，压力项采用PRESTO方法离散，压力、密度、速度、液相率和能量亚松弛因子分别设为0.1、1、1、0.7和0.9。考虑液相区存在自然对流，在材料设置面板密度项选择Boussinesq假设，同时将相变材料区域初始温度设置为335.15 K，时间步长选择0.01 s。为方便分析，进行以下假设：相变材料性质均匀稳定且各向同性；液态相变材料为牛顿流体，液相区流动为二维、非稳态、层流和不可压缩；密度采用Boussinesq假设，其余物性参数与温度无关；计算过程中壁面温度恒定；忽略方腔左右壁面的散热损失，即左右壁面绝热。物性参数见表1。

图1 石蜡方腔模型Fig.1 Physical model of the paraffin in square cavity

表1 石蜡物性参数Tab.1 Physical parameters of paraffin

1.2 数学模型

连续性方程为：

div(Ut)=0

(1)

式中：Ut为速度矢量。

动量方程为：

(2)

式中：τ为时间；Sv为动量方程的源项；ρ为密度；p为压力；ν为运动黏度;u为x方向上的速度；v为y方向上的速度。

能量方程为：

(3)

H=h+ΔH

(4)

(5)

式中：ΔH为潜热；λ为导热系数；T为温度；cp为比热容；Tref为参考温度；href为参考焓。

引入液相率f为：

(6)

式中：L为相变潜热；Ts为开始发生相变的温度；Tl为完全变为液态时的温度。

动量方程的源项Sv为：

(7)

式中：ρref为参考密度；g为重力加速度；β为体积膨胀系数。

1.3 边界条件

下壁面为恒温加热面，则有

y=0,u=v=0,T=Th

(8)

式中：Th为底部加热温度。

上壁面为恒温冷却面，则有

y=A,u=v=0,T=Tc

(9)

式中：Tc为顶部冷却温度。

左右两壁面为绝热面，则有

(10)

引入瑞利数Ra：

(11)

式中：a为热扩散率；Tm为相变温度。

2 模型验证

分别采用不同网格数进行无关性验证，结果见图2。从图2可以看出，网格数为31 773和40 000时液相率曲线较吻合，两者偏差为1%；而网格数为17 888时液相率曲线与前两者的吻合度略差，最大偏差达到2.4%。因此，选择网格数为31 773。

图2 网格无关性验证Fig.2 Grid independence verification

为验证模型的有效性，将本文模拟的液相率与文献[8]中的结果进行对比，见图3。由图3可以看出，本文与文献[8]的液相率曲线基本吻合。

图3 模型有效性验证Fig.3 Validation of the model

3 结果及分析

3.1 多解现象

初始条件设置成模型高度为5 mm，Ra为3.5×104。图4和图5分别给出了第1次模拟时不同时刻下的固液融化图和速度流线。在0～100 s时刻下方腔内石蜡不断融化，液相区高度增加，石蜡液体受到浮升力和重力的影响，自发形成了9个规则且对称性良好的速度涡。由于液相区高度仍较低，Ra较小，自然对流的影响较弱，传热仍以导热为主，固液相界面呈规则的水平状态。随着时间的延长，液相区高度不断增加，自然对流的影响逐渐增强，110 s时固液相界面不再规则。在360 s时刻下，液相区逐渐扩大，Ra增大，自然对流的影响也增强，流场中最左侧的速度涡开始被挤压并最终消失，流场中速度涡的数量减至8个，直至液相率不再改变，达到稳定阶段。

(a)100 s

(a)100 s

图6和图7分别给出了第2次模拟时不同时刻下的固液融化图和速度流线。从2次模拟结果可以看出，其传热过程基本一致，但第2次模拟时初始阶段流场的速度涡数量发生了变化，增至10个，在120 s时刻下固液相界面变成不规则的波浪形，在560 s时刻下流场不再稳定，最右侧的速度涡消失，在640 s时刻下最左侧的速度涡消失，流场中速度涡的数量减至8个，并达到稳定阶段。从最终稳定阶段的速度涡数量看，2次模拟结果相同；从速度涡的方向看，2次模拟结果的速度流线方向相反，这就造成采用相同模型且在初始条件相同的情况下出现多解现象。

(a)100 s

(a)100 s

图8给出了2次模拟的液相率对比。从图8可以看出，初始阶段是以导热为主导的传热阶段，2条液相率曲线吻合度很高，融化速率基本保持一致。在约110 s时刻下第1次模拟的融化速率开始增大，并早于第2次模拟。这是由于第1次模拟结果比第2次更早进入以对流为主导的传热阶段，该阶段的融化速率高于以导热为主导的传热阶段，说明自然对流对传热起到强化作用。随着时间的增加，第1次模拟的液相率始终高于第2次模拟，直至最终2条液相率曲线基本重合。第1次模拟的液相率达到稳定所用的时间少于第2次模拟。2次模拟的液相率达到稳定所需时间和融化速率均存在差异，说明其结果存在最优解。

图8 2次模拟的液相率对比Fig.8 Comparison of liquid phase rate between two simulation results

3.2 流场的非线性特性分析

图9和图10分别给出了模型高度为9.2 mm时的固液融化图和速度流线。在0～<100 s内固液相界面呈规则的水平状态，流场呈规则对称的排涡状态，在100 s时刻固液相界面变得不再规则，自然对流对融化的影响开始体现，在500 s时刻流场的对称性发生了改变，形成的速度涡处于不稳定状态，因此速度涡之间发生了强非线性相互作用，系统进入自组织阶段，1 200 s时刻流场中由初始阶段的9个涡变为大小不规则的3个涡。这是由于部分速度涡失稳破裂，而另一部分相邻的速度涡发生合并，破裂的涡被合并的大涡吸收，速度涡数量减少，流场不再对称。

(a)88 s

图11给出了模型高度为20 mm时的速度流线。随着时间的增加，液相区开始出现一排规则的速度涡，并逐渐变大。在510 s时刻，速度涡发生变化，速度涡之间相互挤压变形并相互融合，最终达到稳定阶段；初始阶段速度涡为规则的排涡，达到1 700 s时刻变成左上和右下各1个速度涡，左下和右上为连贯的速度涡，达到1 710 s时刻左下和右上各1个速度涡，左上和右下为连贯的速度涡，且该现象交替产生。此最终状态说明融化过程中存在传热流动不稳定性，呈复杂的非线性。

选择监测点位置x=10 mm，y=2.5 mm、4.6mm和10 mm。为研究监测点的非线性特性，采用时间序列图法和相轨迹图法进行分析。图12～图14分别为不同Ra下的时间序列图和相空间轨迹图，其中U、V分别为x方向和y方向上的无量纲速度。从图12可以看出，Ra=3.5×104时流场由初始的暂态转变为流动稳定状态，数值解具有稳定定态解，相空间轨迹表示为一定点。从图13可以看出，Ra=2.19×105时流场由初始的暂态转变为周期性振荡状态，数值解为具有周期性的振荡解，相空间轨迹表示为一封闭的曲线(即极限环)，数值解的轨迹在相空间中经过一段暂态过程后均落在极限环上。从图14可以看出，Ra=2.3×106时流场由初始的暂态转变为最终阶段的振荡状态，这种振荡解是非周期的，且具有随机性，在相空间上表示为杂乱无序的混沌状。这说明改变Ra可以控制流场的最终状态。

(a)88 s

(a)100 s

(a)x方向时间序列图

(a)x方向时间序列图

(a)x方向时间序列图

4 结 论

(1)在相同模型和初始条件下进行数值模拟，结果会出现多解现象。

(2)2次模拟的融化速率和液相率达到稳定所需时间存在差异，说明其结果存在最优解。

(3)Ra=3.5×104时数值解具有稳定定态解，当Ra=2.19×105时具有周期性振荡解，当Ra=2.3×106时具有非周期、随机的振荡解，这说明改变Ra可以控制流场的最终状态。