系统比较CMFD软件：在微通道两相流仿真领域的预报效果
积鼎科技基于VirtualFlow数值计算仿真平台，采用先进数值模拟方法可针对气液两相流的多相动力学及其特性开展数值模拟研究。众所周知，对高密度和粘度比的多相流在许多工程应用中发挥着举足轻重的作用，但对此类流动进行稳健和准确模拟的方法仍然难以捉摸。
近年来，工业界始终在同时推动微流动工程应用部件的性能发展和小型化发展，特别在芯片实验室、生物MEMS和微冷却电子设备等领域。在这些部件的微流动通道中，会发生传热和传质过程，可以通过使用多相流来增加传热和传质的过程。更进一步地深入探索两相流机理特性，如界面拓扑结构和压降等方面，可以进一步提高微流动工程应用部件性能的重要控制参数的合理性。
对应用部件中的微流动过程的理解对微流动部件的设计策略有着直接的影响。计算多相流体动力学（CMFD）是一种研究多相流动现象的数值仿真工具，旨在应用服务于微流体部件的设计过程。与较大通道相比，微通道的流动空间尺寸很小，具有可采用DNS方法进行数值仿真的良好尺度特征。由于控制微流动的主要作用力为毛细管力，所以微流动的流动模式更加稳定，且针对微流动开展的数值计算所需建立的计算域相较于一般的CFD特征空间尺度会大大减小。这使得CMFD可以在保障相同的网格分辨率的情况，仍旧能很好地控制计算成本很高。
近年来，CMFD技术在工业界和学界的应用越来越多，CMFD方法中的自由表面跟踪方法如今已成为商业软件中不可或缺的重要组成部分。我们要认识到，各种数值算法和数值模型目前在这一领域还无法提供物理上100%完全可靠的结果。我们还需要进一步致力于客观描述并尽量消除数值仿真在结果上的不足支出以尽可能获得可靠的模拟结果。积鼎科技正在基于目前已相对确定的流动机理对数值仿真开发策略做更宏大的展望，始终致力于建立完善更稳健和更准确的数值仿真方法并推广，目标超越目前微流体部件设计的经验主义范式。
本文采用两个目前较为流行的CMFD软件：F软件与VirtualFlow研究了微型通道中气液两相流的特性。这两款软件具有不同的自由界面追踪（Interface Tracking, IT）能力：F软件采用Volume of Fluid （VOF）方法，而VirtualFlow采用Level Set (LS)方法。VOF和LS是应用最广泛的两种自由界面追踪(IT)方法，目前已成功广泛应用于多相流动的数值模拟领域。Taha等采用F软件中的VOF方法研究了上升的Taylor气泡，并报告了气相中的强再循环模式。Akbar、Ghiaasiaan以及Quian、Lawal也用F软件中的VOF方法研究了水平通道中的两相流。他们的计算工况中，气泡的流动方向与周围液体的运动平行，气泡与周围液体处于并流状态，数值结果显示气泡内未显示再循环流场。Lakehal团队以及Fukagata分别采用VirtualFlow软件以及其他搭载Level Set方法的CFD软件进行了类似的数值仿真研究，在他们的计算结果中都明确发现了由于气泡周围的流动而显著增强的传热。
VOF方法与LS方法作为目前最为流行的两种多相流自由表面捕捉方法，二者应当在充分收敛的计算条件下得到流动规律与仿真质量相对一致的结果。然而，在文献中，我们发现对于分别采用VOF和LS方法的水平管道中的两相流动的数值仿真，得到的仿真结果之间似乎存在明显的差异。
本文的数值模拟使用了两种不同的CMFD软件，分别为某公司的F软件与积鼎科技的VirtualFlow。F软件采用VOF方法，而VirtualFlow采用LS方法。VOF和Level Set作为被广泛使用的两种自由表面追踪方法，其各自有继承了一些有据可查的、积极的和消极的自身的优劣之处，如图1所示。本文的重点不是通过开发修正这些劣势，而是在评估VOF和LS在数值计算过程中的表现时需要合理地重视并尽量规避这些弊端。F软件与VirtualFlow搭载的都是基于有限体积法求解多流体Navier-Stokes方程的压力基求解器。当然，使用的离散化方案和时间积分方法是不同的。表1列出了模拟中应用的方法，VOF（c）中界面的重建是用二阶CISAM格式进行求解的，而LS函数（φ）是用QUICK线性迎风格式进行求解的。
图1:VOF和LS方法各自的优劣之处
本文研究了两种不同的两相流拓扑结构，气泡流和段塞流。气泡流计算结果在图3中展示出来了，根据直观印象，两款软件的计算结果似乎生成了相同的流拓扑结构(见图3(a)-(b))，但是气泡的大小和破裂脱落频率可以直观观察到明显的差异。F软件的VOF方法的模拟结果相对于VirtualFlow的LS显示在较低的脱落频率下产生较大的气泡。进一步观察气泡内和周围流场的流动表明，在VirtualFlow的模拟结果中，气泡内存在再循环流场。该现象究竟属于物理解还是数值上的非物理解还有待通过与实验结果的比较来证明。
（a）F软件计算结果在T=0.0076s时的密度等值线
（b）VirtualFlow计算结果在T=0.0236s时的密度等值线
图3：F软件和VirtualFlow结果的比较气泡流，三分之一的管道被抽出，0−10D。
图3的c图与d图中黑色线为气体积分数c=0.5或LS函数φ=0的等值线，视作气液两相交界线，黑色箭头线为流场流线。
对于段塞流情况（图4），F软件无法提供具有能够形成周期性形成段塞的物理规律的数值解决方案。段塞的产生似乎是由数值扰动引起的，可能是由于表面张力处理产生的杂散流动。段塞流和残差较高的结果（10-4，泡状流，见图6）与自入口 射入的气体射流的破裂过程具有相同的性质。气体射流在入口处形成细长的气体颈部，该气体颈部折断并缩回到段塞中。另一方面，VirtualFlow的计算结果中产生了具有典型特征的周期性段塞流模式。气液界面处的波纹是微通道中段塞流的典型特征，由气泡前后曲率的差异导致的表面张力的不平衡造成。这些特性可以在图中清楚地观察到。这些结果与Chen等人实验结果呈现出的段塞流状态以及Fukagata等和Lakehal等的数值模拟一致。
（a）F软件在T=7.83×10-3s时的密度等值线（破裂后缩成段塞）
（b）VirtualFlow在T=0.0155s时的密度等值线
（c）Chen等关于段塞流型的实验结果
图4：对比F软件和VirtualFlow关于段塞流计算结果(展示1/3管道:5−15D)
为了研究时间积分算法精度格式对计算结果可能造成的影响，我们使用VirtualFlow进行了两个附加的低阶精度格式算法模拟，即Euler一阶和Runge–Kutta二阶积分算法。由于非定常模拟是在0.1的CFL数限制下运行的，因此没有观察到时间积分方案的强烈影响（图5）。
图5：对比VirtualFlow采用不同的时间积分精度格式得到的速度矢量和流向速度等值线
利用F软件对较高残差(10-4)的气泡流进行数值仿真得到的结果（图6）显示，F软件对界面不稳定性的预测对残差具有很强的敏感性，残差将直接影响气体射流破裂的预测结果。计算结果显示气体射流将在入口处形成一个细长的气体颈，形态类似从水龙头滴下的液体，随后气泡将以较低的频率破裂成蠕动的单个气泡，其拓扑结构类似于段塞流。这里给出的结果与图3的a、b图是相同的，都属于非物理的数值现象。
实际上，当界面运动与外部流动方向垂直时，直接用算术平均值进行处理就是最佳的。但我们发现，剪切传递的缺乏不能仅仅依靠平均的算法来全部解释。Ho等人用相场法模拟了水平管道中的气泡流动。他们的研究结果表明，在界面上应用流体物理性质的算术平均值来处理无法捕捉气相中存在的强烈的再循环流场对界面的影响。
对比F软件的VOF方法，VirtualFlow中的LS方法对界面性质却有不同的处理。在VirtualFlow中，使用两个相邻细胞中心值之间的谐波平均值（基于颜色函数）来计算单元的表面粘度。由于F软件作为商软不对外开放其底层代码，所以目前尚不清楚F软件中是否使用了相同的方法来计算细胞表面粘度。在VirtualFlow中，界面宽度是恒定的，最小为一个单元的尺寸。由于颜色函数的平滑变化，粘度在两个相位值之间平滑变化。但F软件的VOF方法 会导致粘度的急剧变化并且会产生更厚的界面（至少两个单元）。因此，界面单元进行数值计算的流体物理性质不同，这将导致界面的厚度会影响相之间的剪切应力的传递。
我们进一步对泡状流的工况提取了气体入口的截面流动情况，以描述两款软件计算得到的流场和界面剪切应力的差异（图7）。通过比较表明，使用F软件可以获得更平滑的速度场变化。除了粘度的处理外，表面张力项的处理在小规模流动中也很重要。此外，针对密度的一致处理对于大密度比的情况变得尤其非常重要，尤其是对于压力校正方程。零发散条件在VirtualFlow中更容易获得（每个时间步长的相对收敛条件为10-6）。F软件中的残差将更难控制，尤其是对于段塞流的情况，即使我们不断调整松弛因子也很难将数值残差调整到比较理想的结果。
（a）两款软件计算得到T=7.8×10-4时入口附近轴向速度对比，黑色线为自由界面
（b）两款软件计算结果在T=7.8×10-4时入口速度在界面处有不同的速度梯度
图7：两款软件在偏转流场和速度剖面方面计算结果的对比
本文对两种CMFD软件F和VirtualFlow进行了两相流：空泡流与段塞流的流型进行了系统比较。数值验证结果表明F软件的VOF方法在预测段塞流的准确度方面存在不足，而VirtualFlow总体表现更好。根据Chen等人的实验，VirtualFlow的LS方法更准确地捕捉到了流动的非稳态模式，从而在定性上同时产生了捕捉到气泡和段塞流的流动物理状态。
数值计算成本方面，F软件和VirtualFlow在规则网格上的计算时间基本相当。但F软件需要更高的近壁网格分辨率，以避免数值干涸，这意味着在实际工程计算中，F软件需要更高的计算成本。
Nourgraliev等人的研究工作表明，当自由表面流动由剪切力主导时，几乎所有的扩散界面模型都无法 正确预测所有的界面不稳定模式的增长过程。本文通过实际算例展示是否能捕捉到剪切力的影响对主要界面不稳定性的增长的预测至关重要，因为其决定了气泡和段塞的形成。我们的结果表明，对界面附近粘度的处理、表面张力项的离散化以及对界面上密度跳跃的一致处理对于获得与微尺度两相流物理实验结果对比良好的数值结果至关重要。根据我们的分析，F软件无法生成段塞流的原因之一可能是其难以确保在每个时间步都达到高精度的零发散条件。
众所周知，正确预测多相流型对于准确估计传热传质过程至关重要。尽管目前针对微尺度传热的CMFD软件仍然不能提供100%完全可靠的结果，但幸运的是，我们仍然可以通过PIV等实验技术来解决学界及工业界对物理不稳定流型的实验数据的大量需求，从而可以进一步验证不同的CMFD软件的预测结果的可靠程度的差异。
