【多相流沸腾03】VOF方法实现原理

为了能够简单高效地处理两相流动问题，Hirt和Nichols于1981年开创性地提出了VOF（volume-of-fluid）方法，并引入了流体体积函数C的概念，流体体积函数C等于一个单元内目标流体体积与该单元体积之比。

C=1表示单元内充满目标流体，C=0表示单元内不含目标流体，而0<C<1表明单元内同时含有目标流体和非目标流体，据此便可以确定相界面的具体位置及形状。

连续性方程和动量方程中的密度和粘度则通过对两相所占的比例进行加权平均得到。

显然，与单相流动相比，两相流动中多了一个变量C，那么根据多一个变量就需要多一个方程的原则，还需要列出一个关于C的方程。不过还好，这个方程是比较容易找的。

因为根据质量守恒的原则，网格内C的增加量就等于流入网格的C（如果存在相变，还需要考虑源项，这里先不予讨论），即

这便是所谓的VOF方程。

方程有了，最关键的问题还是如何求解。

如果采用一般的差分格式，也就是像处理连续性方程和动量方程那样—进行求解，一段时间之后会抹平C的间断性。为了始终保持界面的锐利，Youngs提出了PLIC（ Piecewise Linear Interface Calculation）方法，又称为分段线性界面重构，其思路便是使用直线来近似代替相界面，是商业软件Fluent求解VOF方程的重要方法之一。

PLIC方法

如图，对VOF方程在一个单元内进行积分，有

整理一下，有

对时间采用显式差分格式，得到

也就是说，为了获得下一时间步的流体体积函数C，需要获得C通过每一个单元界面的流量，这与C在网格内的分布，即目标流体在网格内的形状和位置有密切关系。

PLIC方法中目标流体分布的4种类型（法向量nx>0, ny>0）

具体而言，目标流体在网格内的形状和位置分布一共有16种不同组合，上图展示了其中4种情况，其余12种情况则由这4种情况经过对称操作得到。

显然，对于上图情况而言，四种类型流过单元界面的目标流体形状各不相同，这也是为什么要区分目标流体分布类型的原因。

为了判断目标流体分布属于以上四种的哪种类型，可以做如下考虑：首先，网格的法向量

网格相界面与x轴的夹角a为

结合这二者信息，便可以获得不同网格的目标流体分布类型，针对正方形网格的分析结果如下图所示。

VOF方法实现流程

综合以上分析思路，我们便可以获得VOF方法的实现流程：

1. 获得流体体积函数C大于0的单元，计算C的梯度，判断目标流体分布属于哪种类型；

2.根据相界面类型，计算通过单元界面的流量，并更新每个单元的流体体积分数C；

3.根据流体体积分数C，对密度和粘度进行加权平均；

4.计算与流体体积分数C相关的受力，如表面张力等；

5.此时获得的连续性方程、动量方程与单相流动情况一致，按照单相流动的求解方案对方程进行求解即可。

以下结果均为VOF方法计算所得的结果，是不是还挺漂亮的O(∩_∩)O。

当然，以上分析仅仅代表了VOF方法中的一种典型处理方法，其它相关处理方法正在持续更新中，请大家给帕坦卡一些时间。。。