重新审视固-液相变数值分析中广泛采用的焓-多孔介质建模方法

《Journal of Energy Storage》：Reassessing widely adopted enthalpy–porosity modeling approaches for the numerical analysis of solid–liquid phase change

　　

**摘要**

准确模拟固-液相变仍然具有挑战性，因为需要同时表示两种相态及其不断演化的界面。广泛应用于商业CFD软件的焓-多孔介质法提供了一种实用的解决方案，但其预测精度在很大程度上取决于糊状区常数（Amush）的处理方式和网格分辨率。本研究利用针对一种商用生物基相变材料（PCM）的最新实验数据，评估了焓-多孔介质模型的性能，并考察了Amush与网格独立性之间的关系。模拟在12 cm × 5 cm、倾角为90°的矩形腔体中进行，分别从10,000到1,000,000个网格单元进行系统性细化，涵盖三个Amush值。结果表明，低Amush值会导致无法获得网格独立解，并产生非物理的熔化工前沿行为；而较高值虽能提高一致性，但会增加计算成本。为了量化局部偏差，本研究开发了一种基于图像的异或（XOR）分析方法，证实Amush = 1×10?产生的熔化工前沿与实验结果最为吻合。在该值下，该方法的性能接近界面追踪技术，但仍具有较高计算成本。对于所研究的基准问题，这些发现表明了网格分辨率与Amush之间存在一种强烈的且往往被忽视的相互作用。由于这两个参数均影响动量汇在相变区域中的离散表示，在焓-多孔介质模拟中，应对其独立性进行验证而非假设。然而，该相互作用量的强度以及合适的Amush值可能取决于材料、构型、流态和数值实现方式。尽管该方法被广泛使用，但这一相互作用在文献中受到的专门关注很少，也极少被系统性地研究。研究结果表明，传统的网格独立性程序可能不足以满足需求，并强调了基于物理的校准、熔化工前沿解析验证以及重新审视当前固-液相变建模实践的重要性。

**引言**

对熔化和凝固的兴趣日益增长，与学术界开展的大量数值模拟密切相关。这些应用包括航天器热管理、储能设备设计、处理器或光伏面板的电子热管理以及太阳能空气加热器中的热储存[1]–[6]。计算流体力学作为核心方法，用于详细研究传热现象以及更具应用性的课题，如储能（TES）设备的设计，正变得越来越普遍。商业计算流体力学（CFD）软件降低了探索熔化/凝固物理过程的初始门槛。ANSYS Fluent和COMSOL等商业软件包已成为解决这些问题的主要CFD工具，尽管它们所依赖的方法论可追溯至20世纪80年代[7]。对商业软件的依赖导致研究领域出现一定程度的松懈，大量持续工作集中于一种在熔化工前沿建模的核心方法中存在某些已知不一致性的模型上，而这发生在更先进的TES换热器设计不断发展的背景下[8]。

20世纪70年代，针对相变问题的数值模型初步尝试使用了焓模型[9]、[10]和有效比热容法[11]。这些方法忽略了液体运动，是上述模型的关键局限性，因为实验表明自然对流在准确预测中起着至关重要的作用。Sparrow等人[12]是最早将Boussinesq模型与固体PCM耦合以考虑自然对流的研究者之一，但许多其他研究者早在之前的研究中就通过解析方法研究了自然对流[13]–[15]。

追踪固-液边界是液-固相变问题的主要挑战之一。早期的尝试利用了变形网格法，基于固体壁面的能量平衡重新构建网格[16]。Voller等人[17]调查了熔化工前沿追踪方法并介绍了经过测试和比较的新方法。将熔化工前沿建模为多孔介质因其相比其他方法在数值上更为简单而最具吸引力。该方法后来在一项独立研究中被充分开发并与测试用例进行了验证，被命名为焓-多孔介质法[18]。在该方法开发过程中还突出了一项额外优势，即其表面物理合理性，这是基于凝固过程中枝晶间流体流动的压力降测量数据加以证明的[19]。没有任何明显的尝试或讨论将该实验研究的方程形式与其所采用的、在复杂性上差异显著的建模方法进行匹配。

焓-多孔介质模型最早的验证之一由Brent等人[20]完成，基于Gau和Viskanta[21]的实验，研究了热墙和冷墙之间镓的熔化。该模型基于熔化工前沿形状与实验进行粗略验证。其中一项发现强调了即使金属PCM具有高导热率，对流的重要性依然不容忽视。此外，在凝固过程中，镓的晶体学效应导致不规则的熔化工前沿，难以预测。

Lacroix和Voller[22]完成了一项研究，比较了当时两种主导方法——变形网格法和焓-多孔介质法。总体而言，两种方法在比较熔化工前沿时表现相似，在某些情况下一种方法可能比另一种更有优势，特别与网格创建和数值计算量相关。该研究还首次突出了焓-多孔介质法中一个持续的研究焦点——所谓的孔隙率系数（现称为糊状区常数）取值的不确定性。由于孔隙率系数选择不当而产生的两种极端情况是：低值时产生非物理流动行为，高值时缺乏收敛性，而很少考虑所选值如何影响结果的物理准确性。尽管存在这一关于焓-多孔介质法的关键未解问题，其获得广泛应用可归因于其在商业软件中实施的便捷性。

Lacroix和Benmadda[23]使用了与Brent等人[20]相似的数值方法，探索了填充正十八烷的矩形容器中翅片数量和长度的影响。采用了依赖于熔化分数线性压力降关系，其中一个常数未在论文中讨论，可能是针对实验数据调谐的。手绘的实验熔化工前沿与模拟熔化工前沿进行比较以验证数值模型。这些熔化工前沿与实验早期阶段吻合得很好。这些初步结果展示了将人为压力降与局部熔化分数相关联的模型的潜力，但也存在一些与设定压力降幅度和材料特性相关的潜在困难。

Pal和Joshi[24]应用的焓-多孔介质模型与伴随的实验数据进行了比较。任选了1×10?作为糊状区常数，未提供选择该值的任何解释或依据，且将同一值应用于三种不同材料。初步比较显示模型与实验结果之间存在一定一致性，尽管随着实验时间推移这种一致性有所减弱。由于模型性能不是该研究的主要焦点，因此对模型性能的总体讨论有限。

Shatikian等人[25]、[26]完成了一项研究，调查了在恒温和热通量边界条件下翅片腔体的优势。该研究采用了集成于ANSYS Fluent中的焓-多孔介质模型，并指出糊状区常数通常为较大值，范围在10?到10?之间，其模拟使用了Amush = 10?。对该选择进行了有限讨论，尽管注意到了一项使用1.6×10?值的测试，该值由Brent等人[20]和Pal和Joshi[24]建议。结果表明差异不显著。值得注意的是，两篇被引用研究[20]、[24]中的PCM与Shatikian等人[25]、[26]研究中的不同，且未提及模型验证。

Shmueli等人[27]在圆柱形腔体中对Rubitherm 27的熔化进行了实验和数值研究。这项工作代表了较早的糊状区常数调谐示例之一，该过程此后在文献中变得普遍（例如[28]）。该方法涉及使用多个Amush值进行模拟并与实验数据进行比较，通常是实验域内随时间变化的总体熔化分数。该过程揭示了糊状区常数与总熔化时间之间的非线性关系。

这一发现与Ebrahimi等人[29]后来提出的论断一致，即较大的糊状区常数可能导致不稳定性。Ebrahimi等人[29]的研究范围更广，建模了两个理想化问题条件：盒内相变和激光点熔化。他们采用了一种独特的实现网格独立性的方法，基于两个糊状区常数情况之间的相对差异，并确定了一个对糊状区常数不敏感的独立网格。据当前作者所知，这是唯一已知考虑了网格尺寸与糊状区常数之间潜在相互依赖关系的研究。

几项研究扩展了对糊状区常数影响的研究，使用了一组常被引用的实验数据进行验证，即Shokouhmand、Kamkari和Bruno的研究[30]–[32]。这些研究建模了包含月桂酸作为PCM的120 mm × 50 mm腔体[33]–[37]。验证主要涉及熔化分数演化的比较，类似于Shmueli等人[27]所做的，以及熔化工前沿形状的评估，如同早期基准研究。发现的能产生最佳一致性的Amush值往往取决于材料的熔化范围、腔体朝向等因素。在大多数情况下，这些验证努力是在网格独立性研究之后进行的，未考虑潜在的网格相关效应。常见结论得到了强化，包括糊状区常数在自然对流主导的区域影响最为显著。

Kheirabadi等人[38]和Kamkari和Amlashi[40]在COMSOL中使用相同问题条件进行的研究，采用了10?到10?之间的糊状区常数。在给定网格密度下，糊状区常数的影响似乎很小；然而，在非水平倾角下与实验比较不同时熔化工前沿形状时，模型精度下降。Kamkari和Amlashi[40]在自主模型中实现了焓-多孔介质代码，使用了5×10?的糊状区常数。其结果显示出与之前研究相似的困难，在超过水平倾角时难以建模熔化工前沿，且熔化工前沿表现出一些非物理特征，如相对光滑的熔化工前沿上出现明显的尖点，但总体而言与实验熔化工前沿的匹配优于之前的研究。同一模型后来被应用于不同的翅片布置，取得了更大的成功，可能归因于更强的对流和更高的熔化表面积[41]。

同一腔体几何中低于水平面的倾角也进行了数值和实验研究[42]、[43]。这些研究发现，当加热壁位于PCM上方时，一旦熔化进展到加热表面最低点几何下方的区域，熔化速率和熔化工前沿形状在数值上变得特别难以捕捉。对可变材料特性的影响给予了特别关注，并表明这些特性在糊状区内的处理对结果有显著影响。还观察到当对等温壁规定不同温度时，糊状区常数需要调谐。

针对不同几何结构和问题条件的研究也调查了糊状区常数。例如，Arena等人[44]在COMSOL中验证了用于储能的翅片管模型并与实验数据进行比较。

他们发现，糊状区常数在放电（熔化）过程中显著影响温度性能，但在充电（凝固）过程中影响微乎其微，这突出了自然对流在熔化过程中的重要性，从而也凸显了糊状区常数对数值结果的重要性。Hong 等人 [45] 利用 ANSYS Fluent 模拟了低重力环境下矩形几何形状中石蜡的熔化，并将其结果与正常重力环境下球形商业相变材料（PCM）的熔化进行了对比验证。Hong 等人 [45] 建议糊状区常数与糊状区内相变材料固相颗粒的平均直径相关，并呼吁在该领域开展进一步研究。然而，这种解释是对焓 - 多孔介质模型原始应用的一种误解，该模型采用多孔介质方法是因为其“表观物理性”，而非代表糊状区的实际物理机制。近期的一系列研究考察了焓 - 多孔介质方法的几个方面，包括假扩散对模拟精度的影响，这解释了数值模拟与 Gau 和 Viskanta [21]、[46] 先前的实验之间的差异；二维简化在熔化前沿分析中的适用性，这支持了当前的文献实践 [47]；以及引入了一种基于界面位置量化熔化前沿误差的新方法 [48]。同一研究人员的后续研究提出了六水合氯化钙在特定熔化条件下的糊状区常数关联式，显示其对格拉晓夫数和斯特凡数具有明显的依赖性 [49]。开展了几项实验，以更好地检查熔化前沿实验中的糊状区，重点关注厚度、温度和微观结构演化，这些取决于材料和具体问题 [50]、[51]。Yang 等人 [51] 提出了孔隙率常数和平均粒径对液相分数依赖关系的新函数。然而，这些函数尚未在其他研究中得到进一步应用，可能是因为常用商业代码的局限性。

大量的数值研究调查了矩形腔体中的几何变化。Hu 等人 [52] 模拟了正十八烷，以研究楔形腔体中长径比的影响，并使用 Gau 和 Viskanta [21] 的镓数据验证了其结果。Ye [53] 和 Duan 等人 [54] 的类似研究考察了填充石蜡的简单矩形容器的纵横比，发现随着纵横比的增加，熔化时间减少。其他研究人员关注圆柱形几何形状，包括偏心同心管 [55]、用于壳管式换热器的垂直相变材料管 [56]，以及带有铜泡沫插入物的圆柱配置 [57]。Azad 等人 [58]、[59] 的研究也调查了两种不同材料从恒温圆柱的对流起始情况。总体而言，焓 - 多孔介质模型已应用于广泛的几何形状和应用领域，尽管关于适当的验证程序和参数选择仍存在重大未决问题。

在本研究中，使用最近的实验数据 [60] 严格评估了焓 - 多孔介质方法在模拟商业相变材料固 - 液相变中的性能。该研究超越了传统的参数调整，通过系统检查网格无关性和糊状区常数选择的基础，将定量的熔化分数趋势和空间分辨的熔化前沿几何形状直接与实验观测结果进行比较。基于当前作者先前关于热物性定义变化（例如在布西涅斯克近似和温度依赖密度模型之间的选择）可能会显著改变预测熔化行为的工作 [43]，本研究解决了一个额外的且很大程度上未被检查的数值不确定性来源：网格分辨率与糊状区常数之间的内在耦合。尽管焓 - 多孔介质方法在商业计算流体力学软件中广泛实施，并常规接受网格无关性研究，但这两个参数之间潜在的相关性很少被系统地隔离和评估。

该分析考虑了一个倾斜 90° 的 12cm×5cm 矩形腔体，这是先前研究中经常使用的基准配置，以确保与现有文献具有直接可比性。通过明确量化常规的网格细化程序和 ???????????? 选择实践如何影响与实验熔化前沿演化的一致性，这项工作揭示了仅通过整体熔化分数比较所不能明显察觉的现有验证方法的局限性。因此，研究结果为对所调查基准的常见建模假设提供了关键评估，并证明网格分辨率和 ??mush 不应自动被视为独立的数值选择。尽管此处观察到的定量关系无法从单一配置推广，但离散动量汇对网格分辨液相分数场的依赖性提供了一种基于公式的理由，以在其他焓 - 多孔介质模拟中测试这种相互作用。

章节片段
物理模型
图 1 显示了本模型使用的计算域，该模型基于作者之前的实验研究 [60]，并扩展了 Shokouhmand 和 Kamkari [30] 的早期工作。该模型代表一个高度为 ??=120 mm、宽度为 ??=50 mm 的二维腔体，相对于水平方向（＋x 轴）倾斜 90°。为了复制实验条件，在右壁上施加了温度为 57 的等温边界条件。

结果与讨论
二维模型使用 Fidelity Pointwise 2022.1 中的四边形单元结构化网格进行离散化。数值研究的第一步通常是网格敏感性分析。为了启动此过程，使用了 ????????????=1×1?05 的初始估计值（ANSYS Fluent 提供的默认值），并模拟了从 10,000 到 1,000,000 个单元格不等的几种网格分辨率。得出的熔化分数如图 3 所示。当网格细于 50,000 个单元格时

结论
本研究评估了焓 - 多孔介质方法对于包含商业相变材料的垂直矩形腔体的性能。对于此基准和数值实现，结果表明在所研究的较低 ??mush 值下，网格分辨率与糊状区常数之间存在强烈的相互作用。在粗分辨率下获得的表观网格无关解不可靠，因为进一步的细化显著改变了结果。在 ??mush=1×1?09 时，相互作用

CRediT 作者贡献声明
Casey J. Troxler：撰写 - 原稿，方法学，调查，形式分析，概念构思。Nicholas C. Lopes：撰写 - 原稿，方法学，调查。Andrew J. Heiles：方法学，调查。Sandra K.S. Boetcher：撰写 - 审阅与编辑，监督，方法学，调查，形式分析，概念构思。

利益冲突声明
作者声明他们没有已知的相互竞争的财务利益或个人关系，这些利益或关系可能影响本报告中报道的工作。

致谢
C.J.T. 得到了 ASHRAE 资助奖的支持。作者感谢使用安布里 - 里德尔航空大学的 Vega 高性能计算集群提供本研究所需的计算资源。

Casey J. Troxler|Nicholas C. Lopes|Andrew J. Heiles|Sandra K.S. Boetcher
机械工程系，安布里 - 里德尔航空大学，佛罗里达州代托纳比奇，美国