基于LBM的泡沫金属与翅片相变储能系统性能对比分析

张金亚，周文博，程紫漪漪

（中国石油大学（北京），北京 102249）

近年来，可再生能源的可持续性、清洁环保的优点使其受到广泛关注。但由可再生能源的自然特性，在使用的过程中具有间歇性和随机性的缺点，需要与储能系统结合应用[1]。目前储能方式主要分为显热储能、潜热储能和化学储能。其中，潜热储能具有储能时温度恒定的优点，是一种有较大潜力的储能方式[2]。相变材料(phase change materials，PCM)是潜热储能中用于储存能量的介质。然而，工程所用相变材料的热导率相对较低，限制了储能效率[3]。因此，提高相变材料热导率是研究关注的重点。目前，常采用提高潜热储能效率的方法有添加翅片、采用多孔介质、添加纳米材料等[4-6]。

由于制造简单、结构成本低、传热面积大的特点，采用布置内翅片的方法可以显著提高潜热储存单元的有效热导率和传热强度，从而有效提升相变材料的储能效率[7]。由于泡沫金属的高比表面积和高导热性，将相变材料嵌入金属泡沫中制备成金属泡沫/PCM 复合材料同样也是提升相变材料储能效率的有效手段。在金属泡沫/PCM 复合材料中热传递时，热量会沿着金属泡沫基质传递到整个区域，金属泡沫孔隙中的每个PCM 单元都可以被金属泡沫基质加热，从而显著改善PCM的热传导性能[8]。

近年来，学者分别对添加翅片、添加泡沫金属的相变储能系统的传热性能进行了数值模拟和实验研究，Tian等[9]在矩形空腔中插入不同材料的翅片，相变材料的熔化时间均大幅度降低，与径向翅片相比，诸如雪花型翅片[10]、偏心分型翅片[11]能进一步提高相变储能系统的储能性能。Mancin等[12]探讨了3种不同类型的泡沫铜对石蜡-金属泡沫复合材料熔化的影响，结果表明，金属泡沫的存在可以增大相变材料熔化过程中固液相界面的不均匀性，进而改善相变储能系统的性能。Ding等[13]设计了孔隙率为0.95，孔隙密度(pores per inch，PPI，单位英寸长度上的平均孔数)分别为10、20、30 PPI的泡沫金属复合相变材料，使用格子玻尔兹曼方法(lattice Boltzmann method，LBM)模拟了复合相变材料在倾斜45°的方腔中的熔化过程，结果表明，随孔隙密度的增加，自然对流效果降低，平均液相流速和传热系数降低，复合相变材料的熔化时间增加。但在等体积金属情况下，翅片和泡沫金属对PCM 传热效果的强化对比尚未见报道。为此，本工作使用四参数随机生长法(quartet structure generation set，QSGS)构建了孔隙密度分别为20 PPI、30 PPI的泡沫铜复合相变材料模型，并构建了等金属质量的等长双翅片相变材料模型，使用LBM 分别对添加翅片和添加泡沫金属两种结构进行模拟分析，对比研究了其对相变储能系统性能影响的差异。

1 物理模型

1.1 泡沫金属模型构建

QSGS是构造多孔结构的常用方法，其原理是通过在构造区域内随机分布初始固相生长核，让节点上的每个固相生长核以一定的生长概率向其周围点生长，直到达到预定的孔隙率，获得具有一定孔隙度的多孔结构[14-15]。本工作探究的重点在于孔隙密度对复合结构传热性能的影响，为了让生成的多孔结构能够表达不同的孔隙密度，在空间中的固定位置均匀设置初始固相点，固相点的数量即为孔隙密度对应的孔数。通过改变初始固相点的数量改变孔隙密度，进而构建出具有不同孔隙密度的泡沫金属结构，步骤如下：

（1）在空间中设置初始固相点。设置固定位置的初始固相点均匀分布在模拟区域内(10 mm×10 mm)，固相点的数量即为孔隙密度对应的孔数。

（2）按照一定的生长概率(0.01)，令分布的固相单元向8个方向的相邻点生长，随机生成种子概率为0.1。

（3）重复上述步骤，直到生长相达到设定的孔隙率(0.95)时，停止生长。此时重构模型的几何参数有孔隙率和孔隙密度，即QSGS重构多孔介质模型完成。

QSGS 生成的二维多孔泡沫铜结构如图1所示。

图1 翅片与泡沫金属方腔几何Fig.1 Fin and metal foam geometry diagram

1.2 翅片与泡沫金属几何模型

本工作使用的相变材料为石蜡，翅片与泡沫金属的材料均为氧化铜，热物性参数见表1。物理模型如图1所示，方腔的几何尺寸为10 mm×10 mm，选取孔隙密度分别为20 PPI 和30 PPI 的泡沫铜及双翅片结构为研究对象。为了统一金属耗费材料，泡沫金属结构和添加翅片结构中氧化铜的体积均占复合相变材料的95%，即孔隙率为0.95。由于方腔尺寸较小，过多的翅片会使方腔制作困难，因此在左壁面布置等长等距的双翅片，为了降低翅片对自然对流的影响，长度设置为3.85 mm，翅片厚度为0.65 mm，距离上下壁面2.8 mm。翅片与泡沫金属方腔如图1所示，图中黑色区域表示相变材料的填充区域，白色区域表示泡沫金属[图1(a)、(b)]和翅片[图1(c)]。

表1 氧化铜和石蜡的热物理性质[16-17]Table 1 Thermal physical properties of copper oxide and paraffin[16-17]

2 计算模型

PCM 在熔化过程中，受自然对流的影响，液态相变材料向上移动，造成温度、压力分布不均匀。压力和温度与PCM 密度的变化密切相关，使动量方程呈现高度非线性的特征，求解不稳定性增加。因此，引入了适用于浮力密度变化的Boussinesq 近似法来简化方程[18]，即在低速流动中，忽略压力变化对密度的影响，假设只有温度会影响密度变化，浮升力G的计算见式(1)。

式中，β为热膨胀系数；T为PCM的液相温度，K；Tm为相对温度，即相变结束时的温度，K；g为重力加速度，m/s2。

为了简化计算，做出如下假设：①熔化的相变材料是一种不可压缩的牛顿流体。②熔化的相变材料在方腔中的运动为层流，黏性耗散可以忽略不计。③除了由Boussinesq假设引起的密度变化外，流体的热物理特性保持不变。基于此假设，建立描述相变对流过程的控制方程。

格子玻尔兹曼方法是一种介观尺度的计算流体力学模拟方法，该方法通过建立粒子微团简化动力学模型，使用统计平均法获得流体动力学中的宏观量[19]，广泛应用于固液相变的模拟过程。本工作采用Huang 等[20]提出的LBM 固液相变模型，该模型采用双分布函数，分别用于求解速度场和温度场，可以减少迭代步数，节省计算时间。

2.1 相变速度场求解方法

式(2)为速度分布函数的演化方程。

式中，ei为i方向的离散速度，m/s，由式(3)计算，流体粒子移动方向如图2 所示；F为Boussineq 假设中的浮升力；τf为无量纲松弛时间；为平衡态分布函数，对于D2Q9模型来说，其可以通过式(4)计算；B是一个关于流体分数和无量纲松弛时间的加权函数，可以由式(7)计算。

图2 D2Q9格子在二维坐标系上的离散方向Fig.2 Discrete velocities of D2Q9 lattice in 2D Cartesian direction

式中，c为格子声速，即格子长度于格子时间的比值，m/s。

式中，ωi为i方向上的权重系数，其取值如式(5)；cs为格子声速相关变量，m/s，其计算如式(6)；

式中，L为液相体积分数，当L=1 时，B=0，此时PCM 为纯液态，当L=0 时，B=1，此时PCM 为纯固态，L在0和1之间时，表示糊状区的PCM密度分布函数演化方程。

式中，i为格子中离散速度的方向；为i的反方向；us为固相速度，m/s。

采用Chapman Enskog展开，得到宏观密度和速度的计算公式：

2.2 相变温度场求解方法

移动界面上的速度条件采用Noble 和Torczynski[21]提出的浸入式移动边界格式处理，温度场的演化同样包括碰撞和流动过程，温度分布函数的演化方程如式(11)。

式中，τg为无量纲松弛时间；为平衡分布函数，如式(12)：

焓值H由式(13)计算：

温度由总焓决定，如式(14)：

式中，Ts为固相线温度；Hs为相变材料处于固相线温度对应的焓值；Tl为液相线温度；Hl为相变材料处于液相线温度时对应的焓值。

2.3 计算几何模型及边界条件

使用Visual Studio 2013 编制计算程序，方腔为单面恒壁温的强制对流系统。熔化过程中，PCM 的初始温度为303 K，左壁温度恒定为343 K。凝固过程中，PCM 的初始温度为343 K，左壁温度恒定为303 K，其余壁面均为绝热状态，无滑移速度，初始温度与右壁面相等。固液边界采用半步长反弹格式(Zou-He边界)。计算的几何模型及边界条件如图3所示，时间步长为1.0 s。

图3 边界条件Fig.3 Boundary conditions

3 模型验证

3.1 多孔金属相变模拟验证

为了验证构建的PCM-金属泡沫铜系统能否有效预测多孔介质内部的相变过程，基于LBM 模型构建了与文献[22]实验具有相同几何参数(孔隙密度、孔隙率、材料性质和边界条件)的储热系统，模拟了PCM-金属泡沫铜系统中PCM 的熔化过程，记录熔化过程中测点(图4中TC-1位置处)的温度变化，并将结果与Chen的实验结果进行了比较。

图4 泡沫铜石蜡结构及测点位置[22]Fig.4 Schematic diagram of the structure and measuring point location of foam copper paraffin[22]

模拟结果与实验的对比如图5所示，实验与数值模拟的最大误差为6.9%，因此本工作构建的模型能够较好地模拟石蜡在泡沫铜中的相变过程。

图5 泡沫铜-石蜡复合材料测点温度变化曲线与文献对比Fig.5 Temperature change curve of the measuring point of copper foam-paraffin composite material is compared with the literature

3.2 网格无关性验证

选取30 PPI 的金属泡沫模型，模拟石蜡在泡沫铜中的熔化相变过程。设置网格数量分别为40×40、60×60、80×80、100×100，记录液相分数体积随时间的变化率，结果如图6 所示。可以发现3 条曲线基本贴合，熔化速率基本一致，熔化总时间相同，选择网格数量为80×80的模型进行后续模拟计算。

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

4 结果与讨论

4.1 熔化过程分析

图7 展示了熔化时间为1 min、4 min、7 min时的液相率云图。从图中可以看到，在1 min 时，3 种结构的相界面基本都平行于左壁面，且双翅片结构中PCM 的熔化速度最快，这是由于高热导率的金属翅片集中于左壁面热源处，在以热传导为主的熔化初期，高温翅片与PCM 有更多的接触面积。

图7 液相率分布瞬态Fig.7 Liquid fraction distribution transient diagram

在4 min 的时候，在热传导和自然对流的共同影响下，无论是添加翅片还是添加泡沫金属，都出现了不均匀相界面，且方腔顶部的PCM 熔化速度大于底部。添加泡沫金属结构下，方腔顶部的PCM 熔化最快，这是由于添加泡沫金属，会抑制自然对流的发展。而双翅片结构中，在方腔顶部没有任何阻挡物，同时翅片的存在增加了初始熔融区域，促进了自然对流的发展，加速了熔化过程。此外，还可以观察到，30 PPI泡沫金属结构，20 PPI泡沫金属结构底部完全熔化的PCM 面积接近，但由于30 PPI 泡沫金属结构孔隙密度更高，高温金属与PCM 的接触面积更大，所以30 PPI 泡沫金属结构的糊状区面积更大。在7 min 时，双翅片结构下熔化程度高于20 PPI 泡沫金属结构，30 PPI 泡沫金属结构。此时，双翅片结构的上半部分已经完全熔化完毕。20 PPI、30 PPI的固液相界面比之前更不均匀。

努赛尔数(Nusselt number，Nu)是表示对流换热强烈程度的无量纲物理参数，其表示流体层流底层的导热阻力与对流传热阻力的比值。通过计算Nu可以分析3 种结构对PCM 材料熔化过程的具体影响。Nu的计算如式(15)。

式中，Tc为壁温；Hi为方腔垂直方向的高度。

从图8中可以看出，3种结构下Nu都随着时间的增长降低。这是因为在熔化开始时，PCM 和壁面的温差大，靠近壁面处对流换热的阻力小，强度高。随着时间的推进，靠近壁面处的PCM 温度快速增加，PCM与壁面的温差降低，因此Nu迅速降低。而当时间继续推移，更多的PCM 熔化成了液体，此时，PCM 与壁面的温差接近，Nu缓慢降低。

图8 Nu随时间变化曲线Fig.8 The curve of Nu changing with time

熔化开始时，30 PPI 泡沫金属结构的Nu为1.74，远高于另外两种结构。在熔化过程的前3分钟，翅片和泡沫金属的Nu均发生骤降，30 PPI 泡沫金属结构的Nu要高于其它两种结构，这是因为，在换热的前期，泡沫金属结构中的固液相界面较为均匀地分布在壁面附近，而在双翅片结构中，不均匀分的相界面会沿着壁面和翅片分布，翅片的存在会影响PCM 的对流换热强度。随后3 种结构的Nu均缓慢降低，且双翅片结构的Nu最高，20 PPI 泡沫金属结构的Nu最低。这是因为随着熔化的进行，固液相界面移动到方腔中部，自然对流带来的影响高于热传导，而金属泡沫结构中均匀分布的氧化铜会影响对流换热的强度，即泡沫铜结构的Nu都要低于翅片结构，泡沫金属结构的总体换热性能也要低于翅片结构。

图9 是3 种结构在4 min、7 min 时的PCM 速度云图，可以发现，泡沫金属结构中熔化PCM 的最大速度低于双翅片结构。在熔化的全程，双翅片结构下最高速度可达到0.026 m/s，而金属泡沫下最高速为0.017 m/s。4 min时，在方腔顶部，双翅片结构下熔化PCM的速度在0.018 m/s左右，并且液相PCM 几乎覆盖了整个方腔上壁面，泡沫金属结构中方腔顶部熔化PCM 的速度在0.008 m/s 左右，由于金属泡沫的阻挡，对流扩散作用受到影响。7 min时，双翅片结构下，自然对流充分发展，液相PCM的流动范围覆盖了超过3/4的区域，远离壁面处，翅片周围的熔融PCM 流动速度更大，自然对流强度更高，进而在方腔右半部形成了剧烈的涡结构，而在贴近壁面处，速度反而比较小。在泡沫金属结构下，自然对流区域沿着左壁面、上壁面发展，液相PCM的流动范围受到金属泡沫的限制，集中在方腔的左部、顶部。

图9 速度分布云图Fig.9 Velocity distribution program

图10展示了熔化过程中液相PCM体积分数随着时间的变化，从图中可以看出，双翅片结构，30 PPI 泡沫金属结构，20 PPI 泡沫金属结构的整体熔化时间分别为598 s、837 s、725 s。在同样的孔隙率下(使用相同的金属材料)，双翅片结构有最快速的储能效率，而30 PPI 泡沫金属结构的储热性能大于20 PPI 泡沫金属结构，即对于泡沫金属结构来说，孔隙密度越大，对PCM 储热性能的提升效果越低。

图10 液相分数曲线变化Fig.10 Liquid fraction curve variation diagram

因此，双翅片结构对PCM 换热性能的提升作用高于泡沫铜，虽然泡沫铜结构扩展的导热面积更大，但多孔结构会抑制自然对流的发展，从而使得泡沫金属结构下的换热效率低于双翅片结构。如果仅考虑储热情况，不考虑制造方式的区别，在花费同样的金属耗材下，采用翅片结构比采用泡沫金属结构的储能效率更高。

4.2 凝固过程分析

图11 是凝固时间为1 min、2 min 时的液相率分布云图，可以看出在凝固过程中，固液相界面沿着翅片和泡沫金属的形状向右壁面推进。1 min时，由于双翅片结构下翅片分布在冷源壁面附近，热接触面积大，在翅片周围形成了一圈凝固的PCM，此时凝固速度比较快，20 PPI、30 PPI泡沫金属结构的泡沫金属周围同样凝固有固态PCM。在2 min时，双翅片结构中PCM 的凝固体积低于泡沫金属结构。对比1 min 和2 min 的液相率分布云图，可以看到，随着时间的增长，泡沫金属结构的固液相界面较均匀地向右壁面推进，而双翅片结构中，PCM在翅片后的凝固程度高于方腔顶部和底部。

图11 翅片与泡沫金属凝固过程液相率分布瞬态Fig.11 Transient diagram of liquid fraction distribution in the solidification process of fin and foam metal

图12 是3 种结构在1 min 和2 min 的凝固速度分布云图。可以看出，随着凝固时间的增加，凝固的PCM 体积增加，对流换热强度降低，因此在凝固过程中自然对流的影响比较小。在1 min 时，20 PPI泡沫金属结构中的液相PCM运动更加剧烈，最高速度达到0.03 m/s。在2 min 时，三种情况的液体PCM速度都大幅降低，翅片结构下PCM的最高速度在0.012 m/s左右，30 PPI和20 PPI下液相PCM 的最高速度在0.006 m/s 和0.007 m/s 左右，这意味着此时涡流的强度非常小，对流强度较低。在整个凝固过程中，翅片结构的自然对流强度比泡沫金属的自然对流强度高，这与熔化时一致，这是由于翅片结构无遮挡，更加有利于PCM 的对流发展。然而，由于此时三者的自然对流对凝固的影响都比较弱，而且随着时间的影响会越来越弱，因此，即使翅片的自然对流强度更高，总的凝固效率还是要低于泡沫金属。

图12 速度分布Fig.12 Velocity distribution program

图13 是凝固过程中液相分数线随时间的变化曲线，3 种结构下凝固的时间分别为772 s、331 s、264 s。等长双翅片结构的凝固速率远低于泡沫金属结构，在开始时凝固速率较高，随后降低，并以比较缓慢的凝固速度直到相变结束。30 PPI 和20 PPI 泡沫金属结构的凝固速度比较高，30 PPI泡沫金属结构的凝固速度略快一些。

图13 翅片与泡沫金属凝固过程液相分数变化曲线Fig.13 Liquid fraction change curve of fin and foam metal solidification process

在凝固过程中，传热过程以热传导为主，自然对流的影响比较少，尽管翅片和泡沫金属结构具有相同的总金属面积，但由于泡沫金属与PCM 接触的面积更大，且均匀分布在整个储能单元中，使得整个凝固中泡沫金属内部的PCM 都能保持较高的凝固速度。而翅片与PCM 的热接触面积较少，且只集中在方腔左侧，因此，一开始翅片结构的凝固速度比较高，当翅片周围PCM 凝固后，方腔内液相PCM 的凝固速度降低。因此，泡沫金属结构更有利于热传导，也就更有利于液相PCM的凝固。

4.3 储放热总时间对比

图14是3种结构下总储放热时间的对比，等长双翅片结构、30 PPI 泡沫金属结构、20 PPI 泡沫金属结构的总储放热时间分别为1370 s、1168 s、989 s。其中，30 PPI 泡沫金属结构的总储放热时间最短，相比于双翅片结构缩短了27.81%，相比于20 PPI 泡沫金属结构缩短了15.32%。因此孔隙密度更高的30PPI 泡沫金属结构的储放热总效率最高。

图14 相变时间对比Fig.14 Phase transition time comparison diagram

5 结 论

为了对比添加翅片以及添加泡沫金属对相变储能系统换热性能的优化作用，使用QSGS分别构建了等铜体积的双翅片相变储能系统和孔隙密度为20 PPI、30 PPI的泡沫铜复合相变储能系统。采用LBM 对相变材料的储/放热过程进行数值模拟，对比分析了添加翅片以及添加泡沫金属结构对PCM换热性能的影响。结论如下：

（1）熔化过程中，双翅片结构、20 PPI泡沫金属结构、30 PPI 泡沫金属结构的熔化时间分别为598 s、837 s、725 s。双翅片结构的努塞尔数高于泡沫金属结构，即双翅片结构更利于熔化中对流换热的发展，泡沫金属结构中的泡沫铜会抑制对流换热的发展，因此相同孔隙率下采用等长双翅片结构更能提升储热效率。

（2）在凝固过程中，热传导占据主要作用，固液界面分布比较均匀。双翅片结构和泡沫金属结构的涡流强度都明显变低。双翅片结构，20 PPI泡沫金属结构，30 PPI泡沫金属结构的凝固总时间分别为772 s、331 s、264 s。由于泡沫铜的存在大大扩展了PCM 与金属材料的接触面积，因此，相同质量金属材料的前提下，泡沫金属复合相变材料的凝固性能优于翅片结构。

（3）综合考虑储、放热过程，等长双翅片结构，30 PPI 泡沫金属结构，20 PPI 泡沫金属结构的总储放热时间分别为1370 s、1168 s、989 s，30PPI泡沫金属结构的总储放热时间最短，相比于双翅片结构缩短了27.81%，相比于20 PPI 泡沫金属结构缩短了15.32%。因此，综合考虑储放热两个过程，在耗费相同金属材料的条件下，采用泡沫金属结构是更有效的提升储能的手段。