CMFD软件对比：国外商软与VirtualFlow在微通道两相流仿真领域的预报效果

积鼎科技依托 VirtualFlow 数值计算仿真平台，运用先进数值模拟方法，可开展气液两相流多相动力学及其特性的数值模拟研究。在众多工程应用中，高密度和高粘度比的多相流发挥着关键作用，然而，对这类流动进行稳健且准确模拟的方法仍有待完善。本文针对工业界广泛应用的两款多相流数值计算（CMFD）软件 —— 国外商软与积鼎科技的 VirtualFlow，在狭小管道两相流仿真预报方面展开全面对比。两款软件核心差异在于采用不同的常用两相流模拟方法：国外商软采用 VOF 方法，VirtualFlow 则采用 Level Set 方法。对比发现，两款软件预报的流动拓扑结构存在显著不同。对于气泡流，VirtualFlow 的预报结果能够捕捉到气泡内部的再循环流动，而国外商软的 VOF 方法计算结果中未观察到明显的再循环流动。在段塞流场景下，国外商软与 VirtualFlow 关于段塞形成及频率的结果偏差显著。其中，VirtualFlow 的计算结果呈现出周期性的段塞形成，与前人（Chen 2002）的实验结果基本一致；而国外商软基于 VOF 方法的计算结果未能捕捉到周期性的段塞形成过程。

一、引言

近年来，工业界始终在同时推动微流动工程应用部件的性能发展和小型化发展，特别在芯片实验室、生物MEMS和微冷却电子设备等领域。在这些部件的微流动通道中，会发生传热和传质过程，可以通过使用多相流来增加传热和传质的过程。更进一步地深入探索两相流机理特性，如界面拓扑结构和压降等方面，可以进一步提高微流动工程应用部件性能的重要控制参数的合理性。

近年来，CMFD技术在工业界和学界的应用越来越多，CMFD方法中的自由表面跟踪方法如今已成为商业软件中不可或缺的重要组成部分。我们要认识到，各种数值算法和数值模型目前在这一领域还无法提供物理上100%完全可靠的结果。我们还需要进一步致力于客观描述并尽量消除数值仿真在结果上的不足支出以尽可能获得可靠的模拟结果。积鼎科技正在基于目前已相对确定的流动机理对数值仿真开发策略做更宏大的展望，始终致力于建立完善更稳健和更准确的数值仿真方法并推广，目标超越目前微流体部件设计的经验主义范式。

本文采用两个目前较为流行的CMFD软件：国外商软与VirtualFlow研究了微型通道中气液两相流的特性。这两款软件具有不同的自由界面追踪（Interface Tracking, IT）能力：Fluent采用Volume of Fluid (VOF)方法，而VirtualFlow采用Level Set (LS)方法。VOF和LS是应用最广泛的两种自由界面追踪(IT)方法，目前已成功广泛应用于多相流动的数值模拟领域。Taha等采用国外商软中的VOF方法研究了上升的Taylor气泡，并报告了气相中的强再循环模式。Akbar、Ghiaasiaan以及Quian、Lawal也用国外商软中的VOF方法研究了水平通道中的两相流。他们的计算工况中，气泡的流动方向与周围液体的运动平行，气泡与周围液体处于并流状态，数值结果显示气泡内未显示再循环流场。Lakehal团队以及Fukagata分别采用VirtualFlow软件以及其他搭载Level Set方法的CFD软件进行了类似的数值仿真研究，在他们的计算结果中都明确发现了由于气泡周围的流动而显著增强的传热。

VOF方法与LS方法作为目前最为流行的两种多相流自由表面捕捉方法，二者应当在充分收敛的计算条件下得到流动规律与仿真质量相对一致的结果。然而，在文献中，我们发现对于分别采用VOF和LS方法的水平管道中的两相流动的数值仿真，得到的仿真结果之间似乎存在明显的差异。

我们仔细研究了国外商软和VirtualFlow针对微通道两相流的数值仿真结果，旨在回答以下问题：二者的计算结果是否再现了随着气体体积分数的增加而发生的预期流态变化？建立的模型中的哪些差异可以解释数值结果中观察到的偏差？导致入口气体射流破裂的界面不稳定性的产生对自由表面追踪方法数值计算的敏感度如何？并最终通过详细的比较阐述国外商软中的VOF方法用于预测段塞流型的数值方法的不足之处。核心在于其无法预测不断增长的界面不稳定性的形成以及注入气体的最终段塞形成。我们建议或许可以采取缓和策略来改进针对这类问题的建模。

 二、仿真框架

2.1 本文采用的两款CMFD软件的说明

本文的数值模拟使用了两种不同的CMFD软件，分别为ANSYS公司的国外商软与积鼎科技的VirtualFlow。国外商软采用VOF方法，而VirtualFlow采用LS方法。VOF和Level Set作为被广泛使用的两种自由表面追踪方法，其各自有继承了一些有据可查的、积极的和消极的自身的优劣之处，如图1所示。本文的重点不是通过开发修正这些劣势，而是在评估VOF和LS在数值计算过程中的表现时需要合理地重视并尽量规避这些弊端。国外商软与VirtualFlow搭载的都是基于有限体积法求解多流体Navier-Stokes方程的压力基求解器。当然，使用的离散化方案和时间积分方法是不同的。表1列出了模拟中应用的方法，VOF（c）中界面的重建是用二阶CISAM格式进行求解的，而LS函数（φ）是用QUICK线性迎风格式进行求解的。

图1:VOF和LS方法各自的优劣之处

表1：离散化方案和应用的时间积分方法

2.2问题描述

本文针对前人的实验结果的基础上进行了二维轴对称的数值模拟。将空气（ρ_g=1.22kg/m³，µ_g=1.78×10⁻⁵kg/ms）和水（ρ_l=998kg/m³，µ_l=0.001kg/ms）注入一根共流水管中。水管截面的外环圈向管内供水，管轴线上的小同心管向管内供气。假设在t=0时，域中没有完全发展的速度剖面，并且可以忽略入口边缘效应。根据Bretherton推论中的极限，忽略重力效应，Bo<0.842（邦德数定义为重力和表面张力间的比率）。管道的轴向长度为30D，因此完全可以获得完全发展的两相流动。

为了避免国外商软中VOF模型出现的数值干燥（管道壁上液膜的数值破裂，气体与固体壁面直接接触），需要细化近壁网格。使用VirtualFlow的模拟因为没有出现数值干涸，所以可以在等距网格间距的情况下进行。速度和压力的收敛标准均设置为10⁻⁶。边界条件保持与Lakehal等人详细描述的实验相同，初始条件如图2所示。

图2：模拟的初始条件和边界条件。

表2：入口空隙率αi，气体来流速度UG和液体来流速度UL

三、计算技术评估

3.1界面拓扑结构和速度剖面

本文研究了两种不同的两相流拓扑结构，气泡流和段塞流。气泡流计算结果在图3中展示出来了，根据直观印象，两款软件的计算结果似乎生成了相同的流拓扑结构(见图3(a)-(b))，但是气泡的大小和破裂脱落频率可以直观观察到明显的差异。国外商软的VOF方法的模拟结果相对于VirtualFlow的LS显示在较低的脱落频率下产生较大的气泡。进一步观察气泡内和周围流场的流动表明，在VirtualFlow的模拟结果中，气泡内存在再循环流场。该现象究竟属于物理解还是数值上的非物理解还有待通过与实验结果的比较来证明。

（a）国外商软计算结果在T=0.0076s时的密度等值线

（b）VirtualFlow计算结果在T=0.0236s时的密度等值线

（c）国外商软计算得到的T=0.0076s时的压力云图

（d）VirtualFlow计算得到的T=0.0236s时的压力云图

图3：国外商软和VirtualFlow结果的比较气泡流，三分之一的管道被抽出，0−10D。

图3的c图与d图中黑色线为气体积分数c=0.5或LS函数φ=0的等值线，视作气液两相交界线，黑色箭头线为流场流线。

对于段塞流情况（图4），国外商软无法提供具有能够形成周期性形成段塞的物理规律的数值解决方案。段塞的产生似乎是由数值扰动引起的，可能是由于表面张力处理产生的杂散流动。段塞流和残差较高的结果（10^-4，泡状流，见图6）与自入口 射入的气体射流的破裂过程具有相同的性质。气体射流在入口处形成细长的气体颈部，该气体颈部折断并缩回到段塞中。另一方面，VirtualFlow的计算结果中产生了具有典型特征的周期性段塞流模式。气液界面处的波纹是微通道中段塞流的典型特征，由气泡前后曲率的差异导致的表面张力的不平衡造成。这些特性可以在图中清楚地观察到。这些结果与Chen等人实验结果呈现出的段塞流状态以及Fukagata等和Lakehal等的数值模拟一致。

（a）国外商软在T=7.83×10^-3s时的密度等值线（破裂后缩成段塞）

（b）VirtualFlow在T=0.0155s时的密度等值线

（c）Chen等关于段塞流型的实验结果

图4：对比国外商软和VirtualFlow关于段塞流计算结果(展示1/3管道:5−15D)

3.2 VirtualFlow中的高阶时间积分效应

为了研究时间积分算法精度格式对计算结果可能造成的影响，我们使用VirtualFlow进行了两个附加的低阶精度格式算法模拟，即Euler一阶和Runge–Kutta二阶积分算法。由于非定常模拟是在0.1的CFL数限制下运行的，因此没有观察到时间积分方案的强烈影响（图5）。

（a）欧拉一阶

（b）龙格-库塔二阶

（c）龙格–库塔三阶

图5：对比VirtualFlow采用不同的时间积分精度格式得到的速度矢量和流向速度等值线

3.3收敛准则的影响

利用国外商软对较高残差(10^-4)的气泡流进行数值仿真得到的结果（图6）显示，国外商软对界面不稳定性的预测对残差具有很强的敏感性，残差将直接影响气体射流破裂的预测结果。计算结果显示气体射流将在入口处形成一个细长的气体颈，形态类似从水龙头滴下的液体，随后气泡将以较低的频率破裂成蠕动的单个气泡，其拓扑结构类似于段塞流。这里给出的结果与图3的a、b图是相同的，都属于非物理的数值现象。

图6：国外商软计算得到的气液分布(T=0.0076s)证明残差对气泡流拓扑结构具有重要影响

四、对扩散界面单元特性的处理

分析造成上述差异的部分原因可能是两款软件对扩散界面单元特性的处理方式不同。液相和气相之间缺乏剪切应力的传输可能会在流场中引入误差。在界面跟踪方法中，基于界面单元中颜色函数(c_i，j)的分布，在单元(i, j)中计算材料特性(µ, ρ)。VirtualFlow和国外商软都使用算术平均值。

谐波平均在数学上更适用于剪切应力传递描述。

其中剪切应力S_ij给出为，

实际上，当界面运动与外部流动方向垂直时，直接用算术平均值进行处理就是最佳的。但我们发现，剪切传递的缺乏不能仅仅依靠平均的算法来全部解释。Ho等人用相场法模拟了水平管道中的气泡流动。他们的研究结果表明，在界面上应用流体物理性质的算术平均值来处理无法捕捉气相中存在的强烈的再循环流场对界面的影响。

对比国外商软的VOF方法，VirtualFlow中的LS方法对界面性质却有不同的处理。在VirtualFlow中，使用两个相邻细胞中心值之间的谐波平均值（基于颜色函数）来计算单元的表面粘度。由于国外商软作为商软不对外开放其底层代码，所以目前尚不清楚国外商软中是否使用了相同的方法来计算细胞表面粘度。在VirtualFlow中，界面宽度是恒定的，最小为一个单元的尺寸。由于颜色函数的平滑变化，粘度在两个相位值之间平滑变化。但国外商软的VOF方法，会导致粘度的急剧变化并且会产生更厚的界面（至少两个单元）。因此，界面单元进行数值计算的流体物理性质不同，这将导致界面的厚度会影响相之间的剪切应力的传递。

我们进一步对泡状流的工况提取了气体入口的截面流动情况，以描述两款软件计算得到的流场和界面剪切应力的差异（图7）。通过比较表明，使用国外商软可以获得更平滑的速度场变化。除了粘度的处理外，表面张力项的处理在小规模流动中也很重要。此外，针对密度的一致处理对于大密度比的情况变得尤其非常重要，尤其是对于压力校正方程。零发散条件在VirtualFlow中更容易获得（每个时间步长的相对收敛条件为10^-6）。国外商软中的残差将更难控制，尤其是对于段塞流的情况，即使我们不断调整松弛因子也很难将数值残差调整到比较理想的结果。

（a）两款软件计算得到T=7.8×10^-4时入口附近轴向速度对比，黑色线为自由界面

（b）两款软件计算结果在T=7.8×10^-4时入口速度在界面处有不同的速度梯度

图7：两款软件在偏转流场和速度剖面方面计算结果的对比

 五、结论

本文对国外商软与 VirtualFlow 这两款 CMFD 软件在气泡流与段塞流流型的两相流模拟方面展开系统比较。数值验证结果显示，国外商软的 VOF 方法在预测段塞流准确度上存在不足，而 VirtualFlow 整体表现更优。依据 Chen 等人的实验结果，VirtualFlow 的 LS 方法能更准确地捕捉流动的非稳态模式，在定性上可同时准确呈现气泡流和段塞流的流动物理状态。

在数值计算成本方面，国外商软与 VirtualFlow 在规则网格上的计算时间基本相当。但国外商软需更高的近壁网格分辨率以避免数值干涸问题，这意味着在实际工程计算中，国外商软会产生更高的计算成本。

Nourgraliev 等人的研究指出，当自由表面流动由剪切力主导时，几乎所有扩散界面模型都难以准确预测所有界面不稳定模式的增长过程。本文通过实际算例表明，能否捕捉到剪切力对主要界面不稳定性增长的影响至关重要，因为这直接决定了气泡和段塞的形成。分析结果显示，国外商软无法生成段塞流的原因之一，可能是其难以确保在每个时间步都达到高精度的零发散条件。

众所周知，准确预测多相流型对于精确估计传热传质过程意义重大。尽管当前针对微尺度传热的 CMFD 软件仍无法提供 100% 完全可靠的结果，但值得庆幸的是，可借助 PIV 等实验技术，满足学界及工业界对物理不稳定流型实验数据的大量需求，进而进一步验证不同 CMFD 软件预测结果的可靠程度差异。