VOF方法中相变模型综述

摘 要：数值模拟是探究沸腾换热机理的重要方法。本文对使用VOF方法模拟气泡生长的数值方法进行了调研，着重总结了模拟过程中所使用的相变模型，并对相变模型进行了简要的评析，使用其中两个模型模拟了一维Stefan蒸发和冷凝问题，验证了模型准确度。

关键词：VOF方法；相变模型；数值模拟

1 概述

沸腾是指液体内部生成气泡或气相的一种剧烈的汽化过程。沸腾换热则指该过程中的热量传递。沸腾换热由于其优秀的换热能力，广泛地应用于制冷、发电、化工等领域。沸腾过程中伴随着气泡成核、生长、聚并等行为，这些气泡行为体现了两相间的质量、能量、动量传递。对这些气泡行为进行研究，将有助于完善沸腾换热机理，推广沸腾换热的工程应用。

对于两相流模拟，模型主要分为两大类：高相分数模型和界面类模型。前者适用于气泡特别多的流动，不关注气泡界面和气泡形状，着重关注含气率；后者适用于需要捕捉相界面的情况。常用的界面捕捉方法包括Level Set方法和VOF（Volume of Fluid）方法，其中VOF方法由于固有的质量守恒特性，在许多CFD软件中得到应用。由于气泡生长过程中存在相变，因此需要选择合适的传热传质模型结合界面捕捉方法，才能够较为准确地模拟气泡的生长。

从实际需求出发，本文对适用于VOF方法的相变模型进行调研总结，并对不同模型進行了分析比较，以便后续模拟研究中选择合适的模型。

2 控制方程

在相变研究中，研究人员提出了许多相变模型来对气泡生长进行模拟。这些相变封闭模型每个都适合不同的应用，但大多数都采用通用的应用形式，即将相变源项应用于质量、能量和相分数的控制方程［1］。

黏性不可压缩流体的动量方程、质量方程和能量方程如下：

ρu→t+u→·SymbolQC@

u→=-SymbolQC@

p+SymbolQC@

·μSymbolQC@

u→+SymbolQC@

u→T+F→（1）

SymbolQC@

·u→=1ρv-1ρlm·（2）

ρcpTt+u→·SymbolQC@

T=SymbolQC@

·（λSymbolQC@

T）-m·hlv（3）

上式中，ρ表示密度，u→表示速度，p表示压力，μ表示动力黏度，F→包括重力和表面张力，m·表示相变速率，cp表示比热容，T表示温度，λ表示热导率，hlv为液体的潜热。下标v和l分别表示气相和液相。表面张力通过连续表面力（Continuum Surface Force，CSF）模型［2］，作用在相界面处：

Fσ=-σSymbolQC@

·SymbolQC@

αSymbolQC@

αSymbolQC@

α（4）

如前文中所述，需要使用VOF方法对界面进行捕捉。VOF方法使用体积分数α来表征网格单元中的相分布。α=0，表示该网格为纯气相；0<α<1，表示该网格为两相混合，即为相界面所处的网格；α=1，表示该网格为纯液相。求解体积分数α的控制方程如下：

αt+SymbolQC@

·αu→=-m·1ρl（5）

3 相变模型

如上所述，针对相变研究需要选择合适的相变模型来计算相变源项。本节将对不同学者使用VOF方法模拟气泡生长时常用的相变模型进行总结。

3.1 Lee相变模型

根据Lee［3］提出的相变模型，蒸发和冷凝过程传质速率分别表示为：

m·e=reαρlT-TsatTsat蒸发（6）

m·c=rc（1-α）ρvT-TsatTsat冷凝（7）

上式中，re和rc分别表示为蒸发和冷凝的传质强度因子，Tsat为饱和温度。re和rc一般按照经验选取，过大会导致数值收敛问题，取得过小则会导致界面温度与饱和温度的显著偏差。

3.2 Sato相变模型

Sato等人［4］提出了一个简单而直接的相变模型，从温度场计算出的相变速率直接作为体积守恒的源项。相变速率m·定义为：

m·=M·Ai/V（8）

上式中，M·表示界面相变速率（kg/m2s），V表示网格单元的体积。因此对于不包含相界面的网格，其相变速率m·为0。界面相变速率M·定义为：

M·=（λlSymbolQC@

Tl·n→+-λvSymbolQC@

Tv·n→）/hlv（9）

上式中，n→为界面的法向量，由气相指向液相。假设气液相界面温度为饱和温度，用于计算温度梯度。

3.3 Onishi相变模型

Onishi［5］使用温度恢复法，将相变的速率与模拟的时间步长进行关联，得到如下的相变模型：

m·=ρcpTcell-Tsathlv·Δt（10）

上式中Tcell为相界面所在网格单元的平均温度，ρ为相界面所在网格单元的平均密度。Rattner［6］认为应当对式（10）的相变速率进行限制。在一个时间步内，网格单元内蒸发的质量不应超过液相质量，冷凝的质量不应超过气相的质量，对应的蒸发以及冷凝传质速率如下：

m·e，lim=αρlΔt蒸發（11）

m·c，lim=-（1-α）ρvΔt冷凝（12）

类似于库朗数的限制，Rattner认为［6］，为了保证模拟过程的数值稳定性，一个时间步内蒸发产生的气体体积或冷凝产生的液体体积不应超过该网格的体积，对应的蒸发和冷凝速率如下：

m·e，CFL=1Δt1ρv-1ρl-1蒸发（13）

m·c，CFL=-1Δt1ρv-1ρl-1冷凝（14）

因而，模拟相变时，相变速率模型的选取有如下限制：

m·e=minm·，m·e，lim，m·e，CFL蒸发（15）

m·c=maxm·，m·c，lim，m·c，CFL冷凝（16）

对于模拟过程中的时间步长，除了满足CFL条件外，还应满足热扩散稳定条件：

ΔtSymbolcB@

lminΔ2/λρcpeffi（17）

上式中，l为一个自定义的约束，对于二维情况取1/4，对于三维情况取1/6；Δ为网格的最小长度。

4 模型分析验证

上文中所述三个相变模型是目前VOF方法模拟相变时，使用较为广泛的三个模型。Lee模型由于简单而被广泛使用，但它是一种经验模型，模型中的传质强度因子对于不同的工况需要选择不同的值。Onishi模型以及改良后的Rattner模型，采用温度恢复法将相变速率与时间步长关联，但需要对传热速率以及模拟的时间步长进行限制，不然结果将显著失真。Sato模型根据界面两侧的温度梯度差确定相变速率，结果更为准确；但是该方法需要对相界面进行几何重构，得到界面面积、方向、位置等详细的信息，在OpenFOAM中使用该方法较为困难。因此出于后续研究的需求，不考虑使用Sato模型。

目前Lee模型和Rattner模型已经被一些研究人员将其植入开源程序OpenFOAM中。本节将分别使用Lee模型和Rattner模型模拟一维Stefan蒸发问题和冷凝问题，考察两模型的适用性与准确性。

4.1 Stefan蒸发问题

图1为一维Stefan蒸发问题的示意图。初始阶段整个计算域为饱和液体，左壁保持高于饱和温度的恒定温度，受左壁影响，液体开始蒸发，蒸汽层厚度逐渐增加。对于一维Stefan蒸发问题，其解析解为：

δt=2η λvtρvcp，v（18）

上式中η由下式获得：

ηexpη2erfη=cp，vTwall-Tsat πhlv（19）

该解析解作为精确解用于与模拟结果对比，作为相变模型误差分析的参考。

在本研究中，模拟的工质为水，饱和温度为373.15K，热壁温度为373.15K。计算域长度为0.8mm，总模拟时长为1s。模拟中所使用工质的相关物性参数如下表所示。

两种相变模型时，相界面位置随时间的变化以及与精确解的对比。曲线显示数值结果与精确解吻合良好。

4.2 Stefan冷凝问题

图3为一维Stefan冷凝问题的示意图。初始阶段整个计算域为饱和蒸汽，左壁保持低于饱和温度的恒定温度，受左壁影响，蒸汽开始冷凝，液体层厚度逐渐增加。对于一维Stefan冷凝问题，其解析解为：

δt= 2tλlρlcp，l12+hlvcp，lTsat-Twall-1（20）

模拟中左侧为冷壁，温度为363.15K。其余计算参数及物性参数与4.1中相同。

图4分别显示了使用两种相变模型时，相界面位置随时间的变化以及与精确解的对比。结果显示，Rattner模型与精确解吻合较好，Lee模型在初始阶段冷凝速率偏大，导致液层厚度偏大，之后逐渐与精确解吻合。

结语

本文对采用VOF方法模拟气泡生长时所使用的相变模型进行了概述，总结其中十分重要的相变模型，并对模型进行了评析与对比。结果发现，Sato模型较为精确，但其尚未应用到后续研究所需使用的软件OpenFOAM中，Lee模型和改进后的Rattner模型在选择合适的参数时均可得到较为准确的结果。但从实际出发，Lee模型的蒸发冷凝传质强度因子应当对于不同工况甚至不同网格取不同的值，取合适的值较为困难。因此更推荐在后续的研究中选择Rattner模型。

参考文献：

［1］Nabil M，Rattner A S.InterThermalPhaseChangeFoamA framework for twophase flow simulations with thermally driven phase change［J］.Softwarex，2016，5（C）：216226.

［2］Brackbill J U，Kothe D B，Zemach C.A continuum method for modeling surface tension［J］.Journal of Computational Physics，1992，100（2）：335354.

［3］Lee W H.A PRESSURE ITERATION SCHEME FOR TWOPHASE FLOW MODELING［M］.Washington，USA：Hemisphere Publishing，1980.

［4］Sato Y，Nieno B.A New Conservative Phase Change Model for Nucleate Boiling［C］//.20th International Conference Nuclear Energy for New Europe 2011，2011.

［5］Onishi H，Kawamura M，Tada Y，et al.Numerical Analysis on Heat Transfer Characteristics of Looped Minichannel Using PhaseChange VOF Method［C］//.Asme International Conference on Nanochannels，2013.

［6］Rattner A S，Garimella S.Simple Mechanistically Consistent Formulation for VolumeofFluid Based Computations of Condensing Flows［C］//.Asme International Mechanical Engineering Congress & Exposition，2014.

作者簡介：王子威（1998— ），男，汉族，江苏宿迁人，硕士研究生，研究方向为反应堆热工水力。