基于LBM研究骨架对相变材料融化蓄热的影响

张艳勇陈宝明李佳阳

(山东建筑大学 热能工程学院，山东 济南250101)

0 引言

相变材料具有体积小、储能密度大、工作温度变化小等优点，在热能储存和热管理领域具有广阔的应用前景[1-2]。然而，低温相变蓄能材料一直面临着低导热性的主要挑战，阻碍其应用。因此，为了提高相变储能系统的热性能，已经进行了大量研究来增强相变材料的导热性能。Rehman等[3]综述了高导热多孔材料对相变材料强化传热的研究进展；Tao等[4]综述了近十年来潜热蓄热系统中相变材料的发展和性能增强方法。在众多学者的研究中可以看出在相变存储系统中，主要有两种强化传热方法：(1)直接应用高导热插入物[5-7]，如金属翅片、金属泡沫和热管来实现；(2)通过使用微型和纳米添加剂[8-9]，如碳纳米管、碳纳米纤维和金属纳米颗粒，以此间接改善相变材料的热物理性质(如导热性、比热和潜热)。

由于金属泡沫结构具有较高的导热系数和较大的传热比表面积，在相变材料中加入金属泡沫结构是一种很有前途的提高相变材料导热性能和整体传热性能的方法。为了模拟相变泡沫复合材料的相变过程，主要从表征体元REV(Representative Elementary Volume)和孔隙等两种尺度来研究。其中，REV尺度的研究中，通常将复合材料的物性看作整体相同的，孔尺度热流体和几何特征集中在等效流体中，其预测结果显示了整体流体流动和传热特性。泡沫材料和相变材料之间假设局部热平衡或局部热非平衡，分别采用单和双温度模型进行数值模拟。REV尺度虽然能够预测多孔介质的相变过程，但该方法的基本假设失去了微观方面的孔尺度特征。相比之下，基于孔隙尺度的微观研究可以得到多孔介质内部骨架和孔隙之间的详细信息，如局部固液界面形状和传播、金属骨架表面与相变材料之间的温度分布和温差，克服了在REV尺度下多孔介质内部骨架和孔隙结构之间流动和传热方面的不足，越来越受到学者的关注。Hu等[10]基于孔隙尺度，采用数值方法计算了复合相变材料的温度变化和熔化过程，研究了泡沫金属的孔隙率、孔密度等几何参数对复合材料热性能的影响。成骥[11]对单质石蜡相变材料和泡沫铝/石蜡复合相变材料的蓄热性能进行了实验研究，对比了不同工况下的相变响应时间、温度分布和蓄热时间。实验数据表明：相较于单质石蜡，复合相变材料的相变响应时间更短，温度分布也更为均匀。Abishek等[12]通过研究微观形态对金属泡沫石蜡复合材料的影响，发现高孔隙率泡沫金属复合材料有利于提高和控制相变材料的融化速度，可用于热能储存或过程温度控制。陈振乾等[13]及杲东彦等[14]研究了泡沫金属中相变材料的相变熔化过程，发现相变材料中填充泡沫金属，可以有效改善相变材料的温度分布，泡沫金属孔隙率越小，石蜡熔化越快。姚元鹏等[15]通过构建泡沫金属内固液相变传热模型，对方腔蓄热单元中泡沫铜强化石蜡相变蓄热特性进行数值分析，泡沫铜显著改善了石蜡相变的空间均匀性，减小了蓄热区温度梯度，使蓄热速率得到有效提高。

由于多孔介质内部骨架和孔隙之间结构的复杂性，传统的数值模拟方法在计算有关相变材料中添加多孔介质骨架时，往往会遇到内部边界复杂、并行计算效率低等问题，而格子玻尔兹曼方法LBM(Lattice Boltamann Method)在处理多孔介质固液相变时正好弥补了这一不足。杲东彦等[16-17]基于孔隙尺度分析了无量纲参数瑞利数(Ra数)、孔隙率及孔密度等参数对融化相变传热过程的影响，发现相变材料的融化率随着Ra数的增大、泡沫金属孔隙率的减少及其孔密度的增大而增大。Chen 等[18]使用LBM模型模拟了相变材料在泡沫金属中的二维熔化过程，数值计算结果与实验结果定性吻合较好。宋林泉等[19]基于LBM研究了多孔骨架物性参数对固液相变蓄热过程以及糊状区的影响。Tao等[20]采用LBM研究了泡沫/石蜡复合相变材料的蓄热性能发现复合相变材料的传热性能是由导热和自然对流共同决定的。随着孔隙率的降低，复合相变材料蓄热率提高，最大蓄热容量几乎保持不变，但蓄热密度显著降低。He等[21]对LBM在多孔介质中单相和固液相变传热研究中的应用进行了详尽的综述，认为LBM将在多孔介质固液相变传热的研究中发挥越来越重要的作用。文章采用描述糊状区流动特征的多孔介质—多相流复合模型[22-24]及三周期极小曲面方法TPMS(Three Periodic Minimal Surface)生成多孔介质骨架，探究骨架的添加对相变材料融化蓄热的影响。

1 物理模型和数学模型的建立

1.1 物理模型的建立

为了近似描述泡沫材料的内部结构，众多学者对泡沫材料做了许多简化。徐伟强等[25]将泡沫材料看作为规则的立方体结构，简化模型体现了泡沫材料通孔率高、比表面积大等优点，且模型较为规则，网格划分较为简单，计算量较少。但对真实泡沫材料内部结构的描述还不够详细。为了更进一步的描述泡沫材料的内部结构，Sundarram等[26]和Annapragada 等[27]采用面中心法和体中心法描述泡沫材料，但其结构不规则，网格划分较困难，计算量较大且收敛较困难。

采用TPMS生成泡沫材料，不仅能够很精确的控制孔隙的内部结构，而且能产生连通性优异、平滑而连续的孔结构，已应用于许多方面[28-29]。

由TPMS生成的物理模型如图1所示，方腔尺寸长×宽×高＝1×1×1，复合相变孔隙率ε＝1-Vs/Vtotal(Vs为骨架体积，Vtotal为方腔总体积)计算。蓝色为相变材料，白色为固体骨架。相变材料的相变中心无量纲温度Tm＝0.25，相变无量纲温度半径TR＝0.1，即开始融化时无量纲温度为0.15，融化结束完全变成液相时无量纲温度为0.35；多孔介质骨架方腔四周及右壁面绝热，左壁面为高温壁面且无量纲温度Th＝1，初始无量纲温度T0＝0。

图1 物理模型及边界条件示意图

1.2 数学模型的建立

1.2.1 控制方程

基于描述糊状区的两区域模型，有如下假设：(1)液相相变材料流体不可压缩；(2)液相相变材料的流动为层流；(3)固体相变材料和骨架是刚性的；(4)液相相变材料的密度变化符合Boussinesq假设；(5)相变过程中体积变化忽略不计。

孔隙尺度下含糊状区的固液相变的控制方程由式(1)～(3)表示为

式中：u为渗流速度，m/s；t为时间，s；r为液相率(对应于多孔介质孔隙率)，r＝0为固相区，r＝1为液相区，0＜r＜1为糊状区；下标fl和fs分别表示相变材料液相和固相；ρ为相变材料密度，kg/m3；cp为相变材料的热熔，J/K；p为相对压力，Pa；υe为有效运动黏性系数，m2/s；T为温度，K；ktotal为有效导热系数，W/(m·K)；La为相变潜热，J/g；F为外力源项，N。

固体骨架能量方程由式(4)表示为

式中：ks为骨架导热系数，W/(m·K)；下标s为多孔介质骨架。

F为外力源项，由式(5)表示为

式中：υfl为液相流体的动力黏性系数，m2/s；K为多孔介质的渗透率，md；结构函数Fε为多孔介质的形状因子；g为重力加速度，m/s2；β为热膨胀系数，1/K；Tref为参考温度，K。

在糊状区的低含液率区域(r＜rtr)，对应的渗透率、有效导热系数分别由式(6)和(7)表示为

式中：C为糊状区常数；在纯流体区，r＝1，式(1)～(3)转化为纯相变材料固液相变的控制方程；在固相区，r＝0，式(1)～(3)转化为导热的控制方程；在糊状区的低液相率区域，0＜r＜rtr，可看作多孔介质区域，该区域内流动、传热满足式(1)～(3)；在高液相率区域，1＞r＞rtr，液固两相流的表征导热系数、表征运动黏度由式(8)和(9)[30]表示为

式中：ϕ为相变材料固相体积分数，即ϕ＝1-r；υe为有效运动黏性系数，m2/s；υfl为相变材料液相运动黏性系数，m2/s。

相变材料的液相分数由焓法求解，由式(10)表示为

相变材料的焓值En和温度Tfl之间的关系由式(11)表示为

式中：Ens为相变开始时相应温度Tfs的焓值，J/kg；Enl为相变完成时对应温度Tfl相应的焓值，J/kg；TR为相变温度半径，TR＝(Tfl-Tfs)/2，K。

根据相变材料的液相率，r由式(12)表示为

由式(10)～(12)可知，固体骨架传热的纯导热方程(4)中的温度场与液相分数是相互耦合的，可以通过数值迭代求解。为减少控制方程中的变量，探讨相变过程中糊状区对相变过程影响的机理，引入以下无量纲参数有：

其中Fo为傅立叶数；Ste为斯蒂芬数；αfl为相变材料的热扩散率，m2/s；σ为有效热熔比；R为有效导热系数比；Da为达西数。把上述无量纲参数带入到式(1)～(4)中，得到对应的无量纲控制方程分别由式(13)～(16)表示为

式中：σtotal为有效热熔比。

1.2.2 三维相变格子玻尔兹曼方法

在数值模拟中，使用双分布格子玻尔兹曼方法D3Q19 模型模拟在孔隙尺度下自然对流的熔融过程。速度和温度的演化方程及平衡态方程分别由式(17)～(20)表示为[31]

式中：fi(x，t)、gi(x，t)为t时刻x位置的微粒在i方向上的分布函数，简写成fi、gi；δt为格子时间步长；wi为权系数；ei为i方向上的格子速度；cs为格子声速；c为声速。

D3Q19模型如图2所示，模型中的格子声速、权系数和离散速度分别由式(21)和(22)表示为

图2 D3Q19模型示意图

对应源项Fi、Sri分别由式(23)和(24)表示为

速度与温度演化方程中的无量纲松弛因子τf、τT分别由式(25)和(26)表示为

宏观密度、速度和温度由式(19)和(20)求得，分别由式(27)～(29)表示为

F中含有流动速度u，由式(30)～(31)表示为

式中：V为临时速度；c0、c1分别由式(32)和(33)表示为

固液相变LBM由碰撞和迁移两个过程组成，具体步骤如下：

(1)初始化温度和速度分布函数；

(2)在时刻t执行碰撞；

(3)在时刻t执行迁移；

(4)计算出温度，速度，液相率等宏观量；

(5)重复步骤(2)～(4)直到满足收敛判据；

(6)输出相应的宏观物理量。

2 网格无关化及模型验证

在无量纲参数为Ra＝1×105、Pr＝1、Ste＝5、ε＝0.9、R＝10，网格大小分别为60×60×60、80×80×80、100×100×100工况下，计算相变材料液相率随无量纲时间变化曲线，结果如图3所示。可以看出，不同网格下相变材料融化曲线几乎完全重合，因此可以认为选取网格数100×100×100满足网格无关化要求。

图3 不同网格数下液相率V f/V 随Fo数的变化图

为了验证建立LBM模型的正确性，在Pr＝10、Ste＝0.1和Ra数分别为1×104、1×105、1×106工况下，模拟纯相变材料的融化过程。方腔初始温度Ta为0、相变材料融化温度Tm为0、相变温度半径TR为0。方腔左侧高温壁面Th为1、右侧低温壁面Tc为0，其壁面绝热。左壁面平均Nu数由式(34)表示为

Jany等[32]提出的左壁面平均Nu数的关联式由式(35)和(36)表示为

选取y＝0.5的截面，3种工况下分别由式(34)和(35)计算的结果如图4所示。不同工况下与文献[33]的关联式都有很好的吻合度，说明所建立的纯相变三维模型是正确的。

图4 不同Ra下LBM数值解与经典解的对比图

由于模型所研究的问题为添加骨架的相变融化问题，因此必须验证骨架和液相材料之间的耦合问题。将模型简化为方腔内填充骨架的三维自然对流问题，腔体中心加入一个立方体，边长为方腔边长的一半，绝热边界条件。流体Pr＝0.71，Ra数分别取1×103、1×104、1×105、1×106，计算热壁面平均Nu，并与文献[34-35]所得的结果进行对比，结果如图5所示。可以看出计算结果与文献的结果有良好的一致性，侧面说明所建立的模型可以正确的处理骨架与液态相变材料之间的耦合关系。

图5 数值模拟结果与参考值的对比图

3 结果与分析

3.1 无量纲参数Ra 数对纯相变材料融化换热的影响

大量研究表明在纯固液相变中，自然对流对固液相变过程有着很重要的影响，无量纲参数Ra数是体现自然对流强弱的关键参数。故在Ste＝5、Pr＝1、Ra分别为1×104、1×105、1×106时，探讨Ra数对固液相变过程的影响，数值模拟结果如图6所示。红色表示液相相变材料，蓝色代表固相相变材料，介于红蓝之间的为糊状区。在相变初期，3种工况下，介于红蓝之间的糊状区与上下壁面垂直，这主要是因为融化开始时液相相变材料较少，自然对流较弱，热量的传递主要通过相变材料之间的导热。随着时间的推移，相界面逐渐发生弯曲，但同一时刻不同Ra数下弯曲程度不同。随着Ra数的增加，糊状区弯曲越明显，且出现上窄下宽的形状。主要是因为靠近加热壁面处相变材料首先发生熔化，造成相变材料体积膨胀，高温壁面融化的相变材料密度和固相相变材料密度存在差异，相变材料受到热浮升力驱动的影响，液态高温相变材料向上运动到顶部表面，使热量在方腔左上方堆积。Ra数较小时(Ra＝104)，热浮升力引起的自然对流较弱，热量向上方堆积较弱，糊状区规则且弯曲不明显；当Ra数较大时(Ra＝106)，自然对流较强，热量在上方堆积较多，堆积的热量用来融化相变材料，因此糊状区弯曲明显且出现上窄下宽的形状。

不同Ra下相变材料的液相率随无量纲时间Fo数的变化如图7所示，相变材料的融化速率，随着Ra数的增加而逐渐增加，在Ra数为1×104、1×105、1×106工况下，完全融化时间分别为0.38、0.25、0.16。与Ra＝104相比，Ra在105和106时，融化时间分别减少34%、57%。这主要是因为Ra表征自然对流强度的大小，Ra数越大，自然对流越强，融化速率越大。

图6 不同Ra下相场和等温线随Fo数变化图

图7 不同Ra 下液相率V f/V 随Fo数的变化图

3.2 固体骨架的添加对相变材料融化蓄热的影响

在Ra＝1×105、Pr＝1.0、Ste＝5、Tm＝0.25、TR＝0.1、R＝100、ε＝0.95时，相变材料相场图与流线分布随无量纲时间Fo变化如图8所示。其中，红色代表相变材料已经融化，蓝色表示未融化的相变材料固相，介于红色和蓝色之间的部分为熔融态(糊状区)，白色部分为TPMS生成的固体骨架，黑色线型为液相相变材料的流线图。由图8可知，在融化的初始阶段，无论是否添加骨架，固液界面糊状区都与上下壁面近似垂直，说明初始阶段相变传热的方式为热传导。随着融化的进行，比较有无骨架时的相场和流线图，发现糊状区都发生弯曲变厚且出现上窄下宽的形状，添加骨架后糊状区变厚更加明显且变得不规则，主要是因为自然对流热浮升力的存在，使热量向上堆积。高导热率的骨架的添加，热量容易更快的通过骨架传递到相变材料内部，使方腔整体的热导率增强，故糊状区变厚，融化速率更快，方腔内热量传递更加均匀。但高导热率骨架的添加，也会使液态相变材料的速度场的发展受到限制，同一时间下纯石蜡的糊状区弯曲更加明显，流线虽然在方腔内会存在顺时针的环流，但局部的流体流动会受到骨架的阻碍作用，使流线变得不规则且局部形成许多涡流。

在R分别为10、50、100工况下，填充骨架导热系数比β对固液相变融化蓄热的性能如图9所示。由图9(a)可知，骨架的添加能明显提高相变材料的融化速率，缩短融化时间。与纯相变材料融化相比，在模拟的工况下，其融化时间分别缩短了12%、28%、31%。这主要是由于骨架的高导热性能，使热量更容易的沿骨架向方腔内部传播，加快了相变材料的融化速率。高温壁面平均Nu数随无量纲时间Fo数的变化如图9(b)所示，可以看出无论添加骨架与否，融化开始后平均Nu数都迅速下降，然后缓慢下降。对于填充骨架的融化过程，相比于纯相变过程，其平均Nu数都较高，且骨架导热系数越大，相应的平均Nu数也越高。由此可见，填充高导热性能的骨架后，方腔左侧高温壁面的传热能力得到显著增强。在融化末期平均Nu数随无量纲时间Fo的曲线存在明显的转折，转折点正好对应糊状区开始接触方腔右壁面的时刻，这是因为此时相界面糊状区已经发展到右侧绝热壁面，随着相变的持续发展演化，相界面糊状区表面积逐渐减小，糊状区的吸热能力逐渐减弱，导致了融化末期高温壁面平均Nu的降低。

图8 无/有骨架条件下相场和流线分布随Fo数变化图

图9 不同导热系数比R下液相率V f/V 和热壁面Nu avg随Fo数的变化图

复合相变材料孔隙率对相变材料融化蓄热的影响如图10所示。在Ra＝1×105、Pr＝1.0、Ste＝5、Tm＝0.25、TR＝0.1、R＝10以及孔隙率ε分别为0.95、0.9、0.85时，由图10 可知，加入不同孔隙率的固体骨架时，对相变材料融化的影响是不同的。与纯相变材料融化曲线相比，复合相变材料的融化时间分别缩短了13%、18%和24%。分析其原因，主要是因为高导热率骨架的添加，其骨架/相变材料组成的复合相变材料有效导热系数相比于纯相变材料得到大幅度增加，进而强化了相变材料蓄热。此外，复合相变材料孔隙率越低，其骨架部分所占比例越高，有效导热系数越高，所以融化时间减少，融化速率增加，但孔隙率的减小，也会造成相变材料总体蓄热性能的降低，所以要根据实际情况合理选择。

图10 不同孔隙率ε下液相率V f/V 随Fo数的变化图

4 结论

文章基于孔隙尺度对相变材料内加入固体骨架的融化蓄热过程进行了数值模拟，得到的主要结论如下：

(1)采用TPMS方法生成固体骨架，能有效地预测复合相变材料的融化蓄热过程，为基于孔隙尺度预测复合相变材料融化蓄热特性提供一种有效的途径。

(2)对于纯相变材料，随着无量纲参数Ra数的增加，自然对流引起的热浮升力增加，相变材料融化速率增强，完全融化所用时间越短。

(3)固体骨架的添加对相变材料融化蓄热有很大影响，随着骨架导热系数的增加，在骨架和相变材料导热系数比为10、50、100工况下，其融化时间分别缩短了12%、28%、31%。

(4)骨架孔隙率越小，有效导热系数增加，其蓄热速率也大幅度增加，但骨架的添加一定程度上削弱了相变材料的自然对流。