VirtualFlow | 工业应用中的耦合复杂多相流与传热数值计算

多相流以各种形式出现在环境工程和工业工程中，其往往具有相间传热传质的特点，包括自然界或燃烧室中的蒸发分散相、蒸汽轮机中的冷凝、液体喷雾、碳氢化合物输送、可溶性气体吸收等过程。可以看到核反应堆的热工水力学和油气输送可能是两相流存在最多的两个领域。计算热工水力学已经为发生在轻水反应堆瞬态过程中的复杂流动和热/质传输现象提供了越来越复杂且可靠的数值分析。本文旨在提供基于Level Set技术，针对多尺度多相层/湍流动的数值方案。

图1所示的波浪破碎问题涉及不同的长度尺度，作用于不同的时间尺度。波浪在压力和风切应力作用下，从总流中提取动能而形成；波浪最终会破碎分解成小尺度流动，并将动能耗散为热量。波峰翻卷后会形成白浪或微破碎层，作为一个混合两相流区域，该区域的多相流体很难区分。湍流在流动中无处不在，在剪切面的两侧都有湍流流动。大尺度湍流在顶部由气流与波浪夹带的作用形成。核心液体的流动也属于湍流，具有宽带湍流谱。小尺度湍流可能在波峰处演化。由破碎波峰俯冲运动产生的海沫液滴在空气侧湍流的作用下分散并沉积。海底沉积物的输运受到波浪破溃输送湍流的影响。在相间传质方面，海水雾滴通过蒸发导致气-水的通量不守恒，而微小的水波破碎过程则增强了CO2等可溶性气体的吸收。

图1:不同界面长度尺度的波浪破碎问题

自由表面和界面流动是指两相流动问题中涉及两种或两种以上的不相混流体被随时间演变的尖锐界面分开。通常，当界面一侧的流体是在界面上施加切应力的气体时，后者被称为自由表面。界面跟踪/捕获方式能够定位自由表面，其不是通过拉格朗日方法来跟踪界面(例如通过跟踪位于界面上的标记点等)，而是基于欧拉意义上(固定网格)，保持跟踪场的演变来捕获自由界面，例如Level Set函数，体积分数场等。

目前学界主要有两种计算方法来求解捕捉界面流动，可以分为前端跟踪/捕获和接口跟踪/捕获。前端跟踪是指显式参数化界面并及时精确跟踪表面上点的方法，包括Stokes或势流的边界积分方法。其中，在二维或轴对称几何的情况下，界面只是一条曲线，通常通过弧长和时间来参数化。在前端跟踪中，随着界面的演变，标记在拉格朗日意义上被跟踪，尽管底层流场可以在固定的欧拉网格上求解。前端跟踪方法不是本文欲阐述的重点。

接口跟踪方法作为第二种方法，是本文阐述并采用其对多相流动进行求解的重点。接口跟踪方法通过隐式表示捕获接口来处理此类问题，包括Level Set方法（该方法将界面视为在所有空间上定义的函数的水平表面），以及VOF方法，该方法通过跟踪网格中每个计算单元相对于一个流体相的体积分数来捕获界面的位置。体积分数为0或1的单元格不包含界面，而具有分数值的单元格则包含界面。

2.界面跟踪方法  

界面跟踪方法用于预测需要精确界面识别的两相流，例如气泡，液滴或液体射流的破裂等。该方法的关键是使用了具有可变流体性质以及描述表面张力的单流体守恒方程。接下来主要描述VirtualFlow中采用的Level Set方法。在不可压缩流动条件下，用单流体形式表示的耦合流体运动和传热方程采用如下形式：

其中ρ是密度，p是压力，μ是粘度。方程(2)中的右端项表示表面张力及其接触线壁面贡献(即项c)。另外，n表示界面的法向量，κ表示表面曲率，γ表示流体的表面张力系数，δ表示以界面为中心的光滑狄拉克函数。方程(3)中T为温度，Cp为流体的热容，λ为导热系数，Q'''为体积热源。

在VirtualFlow采用的LevelSet方法中，互不掺混的各相之间的界面用连续函数φ表示，φ表示当地到界面的距离，在界面上设为零，在一侧为正，在另一侧为负。通过这种方法，可以识别两种流体，从而将物理界面的位置与连续函数相关联。方程(2)中的流体性质、体力和表面力都局部依赖于Level Set函数，其演化方程(存在相变时)为:

式中，m点为传质速率，其既可以通过界面上的能量跳跃直接确定，也可以使用传热相关性进行建模。如果传质率被迫仅在三线处起作用，则方程(2)中的源项仅应用于包含三线的单元格。如果单位长度的传质率已知，则根据接触线长度和单元面积将其转换为单位面积的传质率。而流体性质，如密度，粘度，热容量和导热系数等则基于Level Set函数进行局部更新，并使用平滑的Heaviside函数在界面上进行平滑:

在实际算例中，Level Set函数在按方程(4)经过一次迭代步后，不再是距离界面的带符号距离。为了恢复其在界面附近的正确分布，将re-distancing方程解为稳态形式:

其中sgn(x)是Signum函数。采用非振荡三阶WENO格式，在(4)的每个迭代步后求解上式。采用三阶龙格-库塔显式时间积分格式求解Navier-Stokes方程和Level Set方程。对流通量采用三阶Quick格式进行离散。采用二阶中心格式差分扩散通量。

3.VirtualFlow多相流软件  

VirtualFlow是一个基于求解多流体Navier-Stokes方程的多物理场、有限体积计算流体力学软件。软件使用结构化网格，但允许将多个块设置在一起。针对多块问题，采用了基于MPI并行算法。网格排列是并置的，因此可以更容易地处理曲线倾斜网格。求解器是基于压力的(投影型)，使用Karki-Patankar技术对低马赫数可压缩流进行校正。可以采用高阶时间推进和对流格式;空间中的三阶单调格式。多相流是用层流和湍流的界面跟踪技术来解决的，在层流和湍流中，流动被认为是一种具有可变物质特性的流体，这些特性根据被流动平流时的颜色函数而变化，从液相中识别出气体流动区域。具体来说，Level-Set和VOF方法(ITM)都可以用于跟踪不断变化的界面。

4.浸入式边界/表面技术  

图2:IST中立方体网格内实体表面(φ_s=0)的表示  

1.电子PCB板的冷却  

工程上往往对风扇对电子电路板PCB的冷却效果感兴趣，特别是对不同特性的风扇放置在不同的位置对PCB的冷却效果的影响。CMFD有助于提供在特定设置下可能发生的情况的数值计算结果。VirtualFlow在该领域旨在帮助重新设计工程操作。得益于浸入式曲面模块，VirtualFlow生成复杂三维算例设置(网格、区域、风扇尺寸、材料属性、边界条件)仅约需要一个小时。

图3：使用IST法对PCB与风扇冷却组件进行网格划分后的数值计算结果  

图3反映了EST模块在简单笛卡尔网格中映射复杂形状组件的能力，并在必要的地方(围绕组件)进行了细化。湍流的k-ε模型结合壁面函数进行稳态解析。网格由125000个单元格组成。如第3.3节所述，只采用高阶方案。该图描述了风扇产生的冷却流对从组件到堆芯流的传热的影响。在这个模块中，可以自动设置风扇的特性，例如指定风扇的数量、轮毂直径、转速和每个风扇的平均流量。薄板元素被人为加厚，以使它们能够由网格表示(至少3个单元应该覆盖内部固体区域)。每个部件都可以指定其热机械性能和内部体积热源。

2.换热器中的两相流换热  

这里讨论了耦合两相流体流动和共轭传热问题，测试案例是处理热管理在热交换器中的应用，特别是在室内空调大型热泵中的使用。本验证算例将一根或一组有限壁厚的管道放置在自由湍流中，代表空气冷却剂。在管道内部，热水通过重力或外部泵送的作用流动。算例比较了管道完全填充或部分填充的两种极端情况。传热将由冷却剂空气向堆芯流动触发，通过三循环机制：在水中对流，通过管道固体壁传导，最后从管道表面向堆芯外对流流动。

图4:使用IST在笛卡尔网格中表示管表面(φ_s=0)  

采用IST方法对装置(单管和组管)进行网格划分，如图4所示。单个管道使用了包含98x20x66个单元格的笛卡尔网格，五个管道组使用了包含138x20x130个单元格的笛卡尔网格。用5个细胞覆盖壁管。如上所述，所有方程均采用Quick格式。Level Set用于跟踪管道内的水面(当它被部分填充时;其中流入孔隙率设置为50%)。在稳态条件下，使用k-ε湍流模型结合壁面函数仅对冷却剂空气进行求解(假设入流湍流强度为5%)在管道内流动的水呈层流状。基于管道外半径（特征长度）的空气流动雷诺数Re=7000。

图5:三维管道内自由表面和传热情况

部分填充单管的模拟结果如图5所示，其显示了流入的冷却剂，管道的自由表面，以及由冷却剂空气热对流而排出的热轮廓。

图6:完全填充管道(a)和部分填充管道(b)流动的温度等高线  

图6显示了管道的二维截面(位于中心)的流动情况，显示了完全填充和部分填充管道中的温度等高线。图6上图显示，由冷却剂外部强制的管道传热在管道周围是对称的，与预期结果相符。事实上，如果流动雷诺数达到一定程度，会引发涡脱落，此时非定常模拟就会显示出一些对称性的丧失。传热的非对称性描绘在图6下图中，其中热量主要向底部被液体占据的区域传递。

3.输油管道中烃类段塞的形成

图7(左)基于BFC的计算多块网格；(右)管道流动的IST网格  

数值模拟算例建立完整三维计算域，如图7所示。管道长6.3m，管径0.14m。如左图所示，采用多块网格策略，用相邻子域覆盖域。这里使用传统的BFC网格和IST网格。网格子域分布在12个CPU上，在Linux PC集群上并行执行MPI。计算网格数量大约30万个。

该仿真基于用于界面跟踪的Level Set方法和用于湍流模拟的LES方法的结合。MILES方法用于模拟亚网格尺度(SGS)扩散效应。入口流速为14m/s，含气率为50%。图8显示两个交叉流位置在一个时间步长的自由液面位移预测，并给出了速度分布情况。

图8:段塞表面破裂和随后的旋涡脱落，颜色表示速度大小  

段塞的形成是由Kelvin-Helmholtz (KH)不稳定性的第一模态增长触发的。在第一模态的尾迹中形成的第二模态，后与第一模态合并最终导致段塞密封。段塞流在距离入口下游6个管径的下游形成。随着时间推移，由液体形成的堵塞会导致更多的液体积聚从而进一步隔离了气体段塞。该机制不断重复发生，导致沿管道形成连续段塞。图9显示了正负涡度水平，表明在段塞破裂前，湍流在段塞的下游已特别活跃。

图9:正涡度和负涡度轮廓(x和z分量)  

4.燃料束子通道分析  

燃料束子通道流动问题属于非常具有挑战性的热工水力学问题。参考Sadatomi等数据，模拟垂直2·3杆束子通道水力平衡两相流的流动特性。实验以常压常温下的水和空气为工质，2·3棒束通道为试验通道，包含6根矩形阵列棒束和2种6个子通道，模拟沸水堆燃料棒束(图10)。在水力平衡流动条件下，获得了各种单相和两相流动下沿各子通道轴线的流量分布数据和压降。装置尺寸如下:杆径:d = 16mm;杆距:p =20mm;节径比:p/d = 1.25;间隙:S11 =S12 = S22 = 4mm;液压直径:Dh =14.3 mm;流过面积:A = 194 mm2。计算域的长度减小到L = 120 mm，管的长度减小到Lt = 20 mm。域的宽度设置为W = 30 mm。通道底部和顶部分别设置流入和流出边界条件。网格由65x65x130个单元格组成。数值模拟在一个8CPU的Linux PC集群上并行运行，耗时24小时。

图10:IST网格和测试通道的横截面几何形状和尺寸  

在预测沸水堆燃料棒束冷却剂的热水力特性时，必须准确地评估子通道之间的流体传递。两相系统中的流体传递由三个独立的部分组成：空隙漂移、导流横流和湍流混合。与前文中的测试算例相同，这里的流动也是多尺度、多特征的，包括流体流动、传热(共轭)、界面流动，同时，流动中的相变过程也至关重要。这里模拟了两种情况:入水速度U_L = 0.1 m/s的层流和入水速度U_L = 1 m/s的弱湍流；流入孔隙率均为50%。图11与图12分别显示了层流于湍流情况下，管与气液界面之间的热分布。

图11:三维和侧视图管道层流中管道温度和气体拓扑结构

图12:湍流工况下管道温度和气体拓扑的三维和侧视图  

由积鼎自主研发的VirtualFlow对多尺度、多流体的层流和湍流气液两相流动有很好的数值仿真效果，其采用最先进的多相流仿真方法。VirtualFlow的亮点在于基于欧拉方法的界面跟踪方法，特别是将Level Set技术用于模拟各种流动。同时本文给出了基于Virtual Flow处理的各种工业实际算例：电子电路板的热管理、管式换热器内三维耦合两相流与共轭换热、圆管内段塞的形成以及子通道流量分析等算例。测试算例评估了应用于两相湍流的界面跟踪方法的预测性能，实际算例的计算结果足以证明对复杂界面两相流问题的数值仿真已是基于VirtualFlow的CMFD解决方案所能解决的。同时，在计算效率方面，本文中的计算算例均可在Linux PC集群上进行数值计算，并能够在24小时内得到计算结果。