基于分级结构骨架相变储热系统强化传热特性

见禹，陈宝明，宫晗语

（1 山东建筑大学热能工程学院，山东 济南 250101；2 山东建筑大学土木工程学院，山东 济南 250101）

如今，随着全球气候变暖、能源短缺等问题逐渐严重，能源供需之间的差距越来越大。为了解决能源供需在时间和空间上的不对称性，提高能源利用率，大力发展能源革命关键支撑的储能技术具有重大战略意义[1]。热能储存包括显热热能储存、潜热热能储存和热化学热能储存，其中利用相变材料（PCM）储能密度大、温度波动小等特点进行潜热储能的效果尤其明显[2]。近年来，固液相变技术的研究与应用越来越广泛，可用于包括建筑节能、太阳能的储存、工业余热回收、锂电池热管理等多个领域[3-6]。石蜡作为一种易于获得、热性能稳定、无毒、无腐蚀性的材料常被用于相变换热[7]，但是相变材料普遍的低热导率始终限制着热量的高效传递[8]。针对这一问题，研究人员已经提出了许多有效方法，包括在相变材料中添加纳米颗粒[9-11]、翅片[12-14]、金属泡沫骨架[15-17]或微胶囊化[18-20]来制备复合相变材料。在上述强化固液相变传热方法中，金属泡沫构建的多孔骨架由于其重量轻、比表面积大和高导热率等优势应用前景广阔[21]。在多孔骨架熔化过程中，通过对距离加热壁面不同位置瞬时温度的研究分析，发现金属泡沫的添加可以使PCM 固相区的传热速率提高5～20 倍，使整体传热速率提高3～10倍[22]。Tian等[23]实验研究发现金属泡沫嵌入PCM 底部加热时，较小的孔径和孔隙率可以更好地改善传热性能。以上研究结果表明，改变多孔骨架不同大小孔隙率和孔径后，金属骨架与相变材料形成的复合相变材料传热效果差异较大。因此，金属骨架对相变传热过程的影响值得进一步分析。在之后的研究中，孔密度[24]和梯度[25]等因素也逐渐纳入研究者们的考虑范围，可以通过孔密度和梯度的加入来改善均匀多孔骨架中由于自然对流引起的熔化不均匀问题。Yang等[26]对采用针翅与金属泡沫混合结构从凝固前沿与凝固温度进行分析，发现最佳传热结构是具有梯度的金属泡沫-针翅混合结构。Ghahremannezhad 等[27]研究了梯度对加热方向热性能的影响，Y方向正梯度金属泡沫采用侧方加热布局与X方向正梯度金属泡沫采用底部加热布局时，对熔化改善效果较为显著。与孔隙率相似，通过对孔密度的梯度分级也可对相变材料熔化过程产生明显影响，Li等[28]通过观察PCM渗入不同孔隙密度金属泡沫的熔化过程发现，顶部高孔密度可缩短相变时间。此外，还可以通过改变方腔的倾斜度[29]、长径比[30]以及泡沫填充率[31]来改善传热效果。结合Wang 等[32]对泡沫金属微观结构的研究发现，有效热导率与泡沫金属的热导率直接相关，这从侧面反映出了传热结构的重要性，骨架分布位置直接决定了复合相变材料的导热性能。改变梯度、针翅结构以及金属微观结构等手段都是对传热结构进行改变从而改变传热效率。

目前对于多孔骨架的研究方向比较多样，无论是梯度泡沫金属结合翅片还是填充率与倾斜度，其基点都是围绕不同孔隙率及孔密度对相变材料熔化效率进行拓展。仅通过均匀填充的形式或是单一的梯度分布对孔隙率和孔密度进行改变时，热量仍然无法较快穿透进相变材料内部。相变界面的传热往往都是上部远快于下部，甚至由于相变材料密度差的存在会使得熔化模型远离加热壁的底部出现熔化死区，难以达到完全熔化。翅片可以有效控制由于密度差产生的自然对流强度，较为粗壮的骨架形态可以保证热量快速穿透，因此对于横纵侧非均匀骨架填充及多孔骨架结合翅片的研究十分有必要。本文参考叶脉结构分布，对不同向的骨架进行非均匀改变，创新性地提出了一种类似于植物叶片茎干的主干与分支两级骨架相结合的骨架模型，通过对此种骨架不同工况的研究总结规律，并在此基础上结合叶脉通道传热结构的启示以及对同心管和翅片的研究，得到了结合横向主干空心骨架形态与翅片优点的新型翅片管骨架，提高了多孔骨架的热传递效率。

1 计算模型和数值方法

1.1 物理模型

仿植物结构是仿生学重要的一部分，通过对抗强风、抗高温环境植物的结构研究，将植物叶脉结构的力学及传热传质特性结合进工程应用中，可以得到更好的物理模型。植物叶脉通过主干叶脉与分支叶脉为整个叶片传输水分和营养，并起到一定的支撑和保护作用。本文基于嵌入均匀金属骨架的复合相变材料，分方向结合叶脉主干与分支的特性，通过在单元模型中加强横向或纵向的骨架结构，构建出一种可用于传热的物理模型，如图1所示。并且此种模型可以结合翅片特点进行研究，采用正六面体穿孔方法构建均匀骨架，并在此基础上对X和Z方向上的骨架粗细分别调整，Y方向作为腔体纵深方向的骨架，暂时控制为定量，得到了多种骨架单元模型，如图2 所示。小模型正六面体边长为10mm，均匀骨架边长为1.222mm，在此基础上进行结构改变或是横纵方向的粗细改变构建了6种骨架结构模型（图3）。骨架整体在Y方向分为两层，每层有6×6 个骨架单元，主要用于研究XZ平面的传热特性，采用50%的粗细变化进行横纵方向的主干、分支骨架区分。多孔骨架模型具有以下特点：①几何结构构建方便，比表面积大而且具有一定的力学支撑特性和传热特性，有应用于工程实际的基础；②对多孔骨架特点进行了表述，符合添加金属泡沫等复合相变材料的传热特点，可以相对真实地反映出复合相变材料的相变过程；③通过对单元骨架模型不同方向参数进行修改，可以得到多种结构的骨架模型，也方便对影响传热的其他因素进行综合研究。采用多重网格对方腔模型进行网格划分，网格单元大部分采用自由四面体，平均单元质量为0.6338。

图1 主干-分支分级结构骨架构建

图2 均匀/横向主干/主干管/径向翅片管多孔骨架单元模型

图3 骨架模型三维示意图

1.2 数学模型

通过焓-多孔介质模型的建立来求解金属骨架方腔的固液相变问题，为了使计算顺利进行需要对相变过程作以下简化：基于有限元方法，流体流动使用层流接口，并将液体石蜡视为不可压缩流体，忽略重力和斯托克斯流惯性项；传热采用流固耦合求解，采用布辛涅斯克近似，忽略流动加热黏性耗散；石蜡和骨架热物性均为常数且石蜡固液相变密度遵循Boussinesq假设。基于上述假设，采用PARDISO法和GMRES 法两种求解方法，结合分离法和迭代法，对动量和能量的控制方程进行了计算。采用隐式时间步长算法和自适应时间步长算法，结合回差公式对步长进行确定，提高了数值鲁棒性。使用GMRES对流体流动和温度进行耦合时，为了避免过早的迭代终止，将误差估计因子设置为20，最大中间迭代次数为100，计算残差容差为0.01。

相变过程的控制方程包括连续性方程、动量方程、能量方程。

连续性方程见式(1)。

式中，ρ为石蜡密度，kg/m³；u、v、w分别为X、Y、Z方向上的速度分量，m/s。

动量方程见式(2)～式(4)。

X方向

式中，p为石蜡的绝对压力，Pa；μ为动力黏度，Pa·s，取4.43×10-3Pa·s；Su、Sv、Sw分别为X、Y、Z轴上施加的源项，具体计算见式(5)～式(8)。

式 中，β为 液 相 率，β= 0为固相，β=1为液相，0 ＜β＜1时为糊状区；α为石蜡的体膨胀系数，K-1，取0.001K-1；T为过程中的温度，K；T0为初始时刻温度，K；Ts为相变起始时刻温度，K；TL为相变终止时刻温度，K。

相变材料传热能量方程见式(9)～式(13)。

式中，cp,l为液体石蜡定压比热容，J/(kg·K)；ΔHl为石蜡潜热比焓，J/kg；hl为石蜡显热比焓，J/kg；hrefl为石蜡基准比焓，J/kg；H为石蜡相变潜热，J/kg。

金属骨架传热能量方程见式(14)。

在石蜡和骨架接触面温度和热量传递表达式见式(15)、式(16)。

式中，TS为金属骨架的温度，K；n为垂直于金属骨架表面方向的坐标，m。

储热量及平均储热速率见式(17)、式(18)。

式中，cp,s为固相石蜡比热容，J/(kg·K)；cp,S为金属骨架比热容，J/(kg·K)；Tc为相变中心温度，K；Tave为相变完成时刻平均温度，K；mPCM为PCM的质量，kg；ms为金属骨架的质量，kg。

平均Nu数见式(19)、式(20)。

式中，Q为复合PCM吸热量，kJ；Aw为热壁面面积，m2；tn为复合PCM熔化时间，s；L为特征长度，m；k为热导率，W/(m·K)。

瞬时液化速率见式(22)。

式中，Vβ为瞬时液化速率；z为时间步长，s。

1.3 边界条件和材料参数

方腔为60mm×20mm×60mm 的长方体，计算域和边界条件如图4 所示，金属多孔骨架孔隙率为0.85；计算模型通过侧方加热，设置左侧为高温边界，初始温度为293.15K，加热壁面温度为328.15K，其余面均为绝热边界。相变材料采用石蜡，金属骨架采用铝合金材质，通过热常数分析仪和差示扫描量热仪（DSC）可测得热导率和石蜡的熔点及潜热等物理参数。不同型号的石蜡与骨架材料参数存在差异，本文综合考虑文献[33]实验参数确定了石蜡和铝合金材料参数，分别见表1和表2。

表1 石蜡物性参数

表2 骨架物性参数

图4 含金属骨架相变腔体物理模型

1.4 模型验证

1.4.1 模型可靠性验证

本文主要研究方腔中石蜡熔化相变传热问题，采用Kamkari 等[34]月桂酸相变熔化的研究结果和Zhuo等[35]温度变化的研究结果对本文的数值模拟程序进行了可靠性验证。通过图5(a)、(b)将液相率和温度计算结果分别与文献中的模拟与实验结果对比，发现液相率随时间的变化趋势一致，而且温度的数值模拟结果吻合较好，模型验证误差均在±2%以内，从而表明了本文数值模拟的可靠性。

图5 模型验证

1.4.2 网格无关化与时间独立性验证

网格与时间步长对模拟结果是很重要的影响因素。如图6和图7所示，分别将自由四面体网格数量划分为470328 个、1206559 个、2355219 个，并对1s、2s和5s三种时间步长进行了验证。从图6和图7中可以看出，不同网格数量与时间步长得到的结果差异非常小，误差均小于1%。基于计算成本的考虑，本文选用时间步长2s、网格数量为1206559个。

图6 网格独立性验证

图7 时间独立性验证

2 结果与分析

2.1 固液相变界面演化

为了对内部不同骨架之间的方腔传热效果进行研究，分别对200s、400s、600s、800s工况时的三维相变熔化界面进行讨论，见图8。工况1 为均匀骨架、工况2为横向主干骨架、工况3为纵向主干骨架、工况4为横向翅片骨架、工况5为横向主干管骨架、工况6为径向翅片管骨架。首先通过工况1～3 对主干-分支骨架和均匀骨架的效果进行对比研究，从相变界面演化情况来看，工况2的传热效果远远优于工况1、工况3，采用主干-分支骨架结构对均匀骨架的传热效果改变明显。在侧方加热情况下，增大横向骨架体积分配可以加强热量向相变材料内部穿透的能力，通过对均匀骨架体积的重新分配可以得到更有针对性传热的横向主干骨架。接下来以横向主干骨架为基础进一步讨论，将加粗后的主干骨架看作圆管，那么就可以结合管壳式强化传热方面的研究，同心套管结构符合叶脉通道内部空心传质的仿生结构，同时结合翅片管起到增强导热的目的。工况4和工况5分别对横向翅片与主干管进行研究，工况4由于添加了横向翅片增大了骨架与相变材料接触面积，在200s 时可以明显观察到相面熔化界面较工况2 要超前，但是在600s 之后，工况2的效果突显，这是因为虽然主干骨架上的翅片在熔化前期可以加快相变速率，但是熔化后期由于对流作用的增强以及后端翅片距离加热端较远，导致横向翅片传热效果不佳。工况2 与工况5相变情况相似，但在进一步扩大主干管径为2.2mm的测试中传热效果差异可达4.8%，主干管可以变相地扩大传热接触面面积，所以保留主干管结构特点。既然横向翅片不利于液相石蜡传热，那么工况6尝试添加径向翅片来避免影响横向流动，在主干管基础上添加径向翅片，工况6在熔化前期保留了翅片快速传热特点，在翅片周围存在更快的相变趋势而且管内流体率先熔化成为一部分热源（图9），因此当熔化进行到后期时依然起到改善传热的作用。

图8 不同骨架结构固液相变界面演化

图9 t=400s时工况6熔化前沿局部放大图

2.2 不同结构对流速分布的影响

相变传热过程是一个复杂的耦合过程，在流固耦合过程中受骨架影响较大，而且导热效果与对流传热效果相互关联相互影响，为了研究相变过程中通过对骨架导热的调整而产生的对自然对流的影响，下面分别选取200s、400s、600s、800s对熔化初中后期流动情况进行表述，绘制流速矢量图并添加体箭头来直观表示不同结构方腔中的速度分布及流动方向。如图10所示，工况1为均匀骨架方腔熔化过程中的速度分布情况，熔化初期靠近加热壁面的石蜡率先熔化，并且在密度差作用下高温液相石蜡向上运动，靠近固相一侧也产生了相对较小的速度，形成了环流运动趋势，随着时间进行，环流逐渐强烈，并且在方腔上部产生了突起，速度交界面倾斜程度逐渐加大，这也表明了自然对流的加入程度，总体熔化趋势呈环状流动伴随着熔化前沿处的逐渐倾斜。在工况2时，由于横向主干加粗，环流产生的速度加快，200s 时就可以看到明显的环流现象，并且速度在横向主干处向相变材料内部延伸，说明横向主干周围的固相石蜡已经超前熔化，液相石蜡快速流动，整个环流不仅在方腔顶部出现突进，在横向主干附近同样产生向前突进现象。工况3向相变材料内部传热能力较为有限，流体流动速度缓慢，虽然同样会形成环流但是剧烈程度大幅度减小，熔化前沿倾斜不明显。工况4、工况5 与工况2相似但是工况4横向主干附近的突进速度较小，原因是横向翅片对自然对流的阻碍。工况6从初期横向骨架附近的速度突进比较明显，并且熔化前沿处速度倾斜程度大，环流能较快速的朝方腔下部运动，800s 时由于整体方腔大部分熔化完成，速度急剧降低。除了总结不同工况的流动情况差异，各工况流动状况也有共通点，在方腔角落，尤其是左下角会产生朝沿对角线方向的一道速度骤降区域，石蜡流动速度越大越清晰，可以理解为快速流动的液相石蜡在此处由于流向转弯而产生了堆积现象，这是采用方腔对流动的一个不利因素。

图10 不同结构骨架不同时刻速度分布矢量图

2.3 不同结构液相率曲线及瞬时液化速率

为了更好地了解整体的相变规律，现通过液相率曲线来进一步研究不同工况复合相变材料熔化情况，如图11所示。工况3完全熔化时间最长，说明严重抑制了整体的传热速率，即纵向主干结构明显不利于热量向内传递，在其余工况中工况2、 4、5、 6 液相率曲线比较陡峭且完全熔化时间均小于工况1，说明以横向主干为基础可以很好地提高相变材料热性能，而且其中工况6最为突出，即径向翅片管骨架，其完全熔化时间为822s，相较于均匀骨架1314s，完全熔化时间缩短了37.4%。

图11 不同工况液相率随时间变化

图12为不同工况复合相变材料瞬时液化速率随时间变化。通过计算不同时刻液化率变化与时间步长的比值，可得到瞬时液化速率。在初始时刻由于加热面的存在，所有工况的瞬时液化速率都很高，但是随着向内部逐渐加热，距离加热面越来越远，石蜡方腔液化速率会越来越慢。将整个液化速率曲线分为两部分研究，即初始的骤降段以及后期相对缓慢下降段。在熔化开始后，加热面影响骤降，骨架开始发挥作用，这时就处于熔化初期骨架主导传热的阶段，此时可以明显看到横向主干类型的骨架液化速率均位于均匀骨架之上，在降低到一定位置后，会出现第一个波谷，并且波谷在图中出现时的位置越高，代表此工况在初期的液化速率越高，工况2、 4、 5、 6波谷分别位于图中A、B、C、D点，均匀工况波谷高度为0.00105，C点瞬时液化速率为均匀工况的1.88倍。第一个波谷的出现代表已经达到了初期熔融石蜡和骨架传热共同作用的瞬时液化速率，虽然熔融石蜡不一定开始大幅度流动，但是已经完全液化的石蜡已经拥有了一定向固相石蜡传热的能力，这一部分传热的加入以及已经被加热的骨架对相变材料内部延伸的导热使得在第一个波谷之后出现了一个波峰。对于中后期近似的线性下降段，速度衰减越快不仅仅表明传热能力下降，也代表了熔化前沿已经推进到方腔后段，所以在同等层次，相对缓慢下降段的坡度越大，整个方腔的换热能力就越强。工况6在液化速率骤降段曲线最稳定，波谷位置最高，且在相对缓慢下降段坡度最大，从整体瞬时液化速率曲线来看优于其余工况。

图12 不同时刻瞬时液化速率随时间变化

2.4 内部温度均匀性及传热特性

通过方腔熔化过程中相变材料温度最高点与最低点的差值求得每个时刻的温差。从图13 可以看到，总体温差变化分三个阶段，第一阶段骤降段，在热壁面作用于方腔时，热量通过骨架快速传导，使得加热壁面与相变材料内部温差迅速缩小；到熔化中期，靠近热壁面一侧围绕金属骨架附近的石蜡大部分已经熔化，环流形成，石蜡不断循环通过潜热蓄热达到相变温度液化的过程，将熔化界面向前推进，此时就处于第二阶段缓冲段，缓冲段占据了整个熔化过程较多的时间，导热与自然对流传热共同推进方腔内的相变材料稳步熔化，传热速度较快，但是大量热量通过潜热储存进液态石蜡，温差效果难以显现。随之出现第二次温差骤降，这代表熔化过程已经进入末期，在方腔角落的最后一部分固相石蜡已经开始液化。通过对比完全熔化时间发现，在达到完全熔化状态时方腔内仍然存在温度差，这是由于液态石蜡也分为高温熔融石蜡和刚液化时的低温熔融石蜡，并不是达到完全熔化状态温度就会完全均匀。在整个熔化过程中，工况6处于温差最低侧，而且在三个阶段切换过程中，工况6曲线趋势最先达到下一个阶段，这说明径向翅片管结构更加合理。在800s处，工况6已经完成了在缓冲区内的传热，此时工况1均匀骨架内部最大温差为16.88K，工况6 内部最大温差为13.96K，工况6温差相较于工况1缩小了20.9%。

图13 不同工况复合相变材料内部最大温差变化

图14 通过Nu随量纲为1 常数Fo的变化反映加热面附近壁面处对流作用强度。从图14 中可以看到，随着石蜡开始熔化，流动作用快速加强，然后随着石蜡环流逐渐向内部移动，壁面处的自然对流持续衰减。在熔化前期，工况2 的自然对流作用最为剧烈，因为熔化前期传热以石蜡和骨架的导热为主，粗壮的横向主干骨架热传导较为迅速，使壁面处熔融石蜡快速形成了环流，由于空心管内容纳了一部分石蜡，并未直接参与自然对流的运动，所以工况5 强度稍有减小。随着时间进行，固液相变界面向内推进，熔融区厚度增加，壁面处的流动阻力增加，Nu迅速减小，此阶段只有工况4 流动更加活跃，说明横向翅片将流动抑制在加热壁面附近，阻碍了环流向内推进。总的来说，横向主干骨架对流动促进作用最明显，添加翅片等其他因素会影响流动效果，而且对同时添加翅片及管结构的工况6 影响更为显著，但是影响流动并不一定传热效果不好，结合图13温差变化图来看，工况6 虽然抑制了对流传热的效果，但是总体的传热均匀性反而最好，这得益于导热和流动传热的共同作用，相变传热过程并不是一个简单的叠加过程，对于单纯强化流动传热或导热都不是最佳选择，针对本文复合相变材料方腔模型，虽然在相变传热中很大一部分是靠对流来实现的，但是适当抑制对流反而能达到更好的传热效果。

图14 不同工况对熔化过程中平均Nu的影响

以相变完成时刻的方腔内部温度为基础，可以计算不同方腔的储热量，并结合各工况的完全熔化时间可以得到各方腔的平均储热速率，见图15。在定孔隙率情况下，每种工况的骨架石蜡体积都相同，但是储热性能差异十分明显。储热量受方腔内部传热影响，在温度均匀性部分也提到过，液相完全并不代表温度均匀，受时间和温度均匀性影响，各工况储热量略有差异，但相差不大。储热速率受完全熔化时间与内部温度梯度的影响较为明显，当温度通过骨架快速将热量带入相变材料内部时，相变材料内部产生较大的温度梯度，使得完全熔化时间缩短，从而体现出储热速率的差异。工况6为强化传热速率的最佳工况，径向翅片管工况内部温度梯度较大，完全熔化时间最短，平均储热速率较均匀骨架可提高61.9%。

图15 不同工况储热性能图

3 结论

本工作将主干-分支、翅片等不同结构骨架加入石蜡方腔之中制备6种工况的复合相变材料，采用有限元方法模拟方腔内的石蜡相变过程，通过对相界面演化、液相率、传热特性等结果分析，得出以下结论。

（1）采用主干-分支骨架对相变传热过程影响显著，在提出的6种工况中，径向翅片管骨架保留了横向主干穿透传热以及翅片扩大接触面的特点，加快了相变界面的移动。

（2）熔融石蜡形成了环流流动，含径向翅片管骨架形成的环流不仅从方腔顶部突进，在横向主干附近同样会产生流速突进现象，速度前沿倾斜程度大，能快速移动至方腔下部。各工况在方腔角落尤其是左下角会产生石蜡堆积，不利于流动。

（3）径向翅片管方腔完全熔化时间最短，在液化速率骤降段曲线最稳定，波谷位置最高且在相对缓慢下降段坡度最大。相较于均匀骨架，波谷高度增加了1.88倍，完全熔化时间缩短了37.4%。

（4）在缓冲区大量热量转化为潜热，温差降低效果不明显，而且达到完全熔化状态时方腔内依然存在温度差。横向主干空心径翅片工况在800s缓冲区传热完成时的温差相较于均匀骨架减小了20.9%。

（5）传热过程不是简单的流动传热与导热相叠加的过程，适当抑制对流反而能提高温度均匀性达到更好的传热效果。受熔化时间与温度均匀性影响，各工况储热量略有差异但相差不大。含径向翅片管骨架的方腔熔化时间最短，与均匀骨架方腔相比，平均储热速率可提高61.9%。