interFoam进行两相流模拟的一些典型参数设置对比&paraview绘制相界面的方法

1 背景

我的个人课题是基于interFoam的一些二次开发,然而我在使用自己的求解器进行模拟的时候往往会因为模拟参数配置的不合适导致这样或那样的一些问题。这些问题本质上是自己对interFoam进行模拟时的一些关键参数不太了解导致的。前几天看到一本资料上关于interFOAM的讲解，其中的拓展练习中包括一些不同模拟参数的对比，借此机会我就使用溃坝算例对一些典型参数对模拟的影响进行了实际的试算分析,也将之前积累的一些心得一同分享出来。

2 算例

二维溃坝算例，来自OpenFOAM7的自带的教程。将模拟参数延长为2s而已。

3 参数对比

3.1 压力边界条件？

固壁处的压力边界条件最常用的是零梯度边界条件，但是当求解器中植入了额外的体积力时（比如两相的表面张力和重力），零梯度边界条件会造成边界处的通量不守恒。OpenFOAM官方手册对fixedFluxPressure是这么描述的：

“At all wall boundaries, the fixedFluxPressure boundary condition is applied to the pressure field, which adjusts the pressure gradient so that the boundary flux matches the velocity boundary condition for solvers that include body forces such as gravity and surface tension.” Manual 2.3.2节

“In the above example, it can be seen that all the wall boundaries use a boundary condition named fixedFluxPressure. This boundary condition is used for pressure in situations where zeroGradient is generally used, but where body forces such as gravity and surface tension are present in the solution equations. The condition adjusts the gradient accordingly.” Manual 5.2.3节

在溃坝算例中，如果我们把fixedFluxPressure边界条件改为零梯度边界条件，那么其会直接由于初始通量不守恒无法开始运算，求解器会要求你使用potentialFoam对速度场进行初始化。然而由于溃坝算例没有进口类和出口类边界条件，因此其实质上无法进行初始化。对于有进口类和出口类边界条件（轴对称和周期性边界条件也算）的算例，此时使用零梯度边界条件是可以进行计算的，但是会造成较大的误差，在壁面附近可以观察到比较明显的误差区域。

3.2 相分数对流项的离散格式？

在默认的算例中，我们的想分数默认离散格式为：

    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss linear;

vanLeer是具有TVD性质的离散格式，一般来说其可以在解变化比较尖锐的区域也有较好的求解精度。在这里我们可以将其替换为一阶迎风格式，两个格式的计算结果对比如下：
upwind格式最大的缺点是其数值扩散过强，这点在结果中也有所体现：右侧图片中相界面外侧低相分数区域的尺寸upwinde格式要更大一些（动图看不清楚的同学可以自己模拟后一帧一帧对比看看）。

3.3 界面压缩因子？

界面压缩因子用来设置interFoam相分数方程中人工对流项的强度（详见李东岳老师的interFoam解析）。如果界面压缩因子c设为1则表示人工对流项完全起作用，对界面进行法向压缩。如果界面压缩因子c设为0，则表示界面压缩的人工对流项完全不起作用。在这里我们分别测试了c=0和c=1的两种情况：
可以看到，人工对流项不存在的时候界面附加出现了严重的数值扩散，界面变得非常“厚”。在坝体右侧区域甚至出现了相分数完全处于0与1之间的情况。

3.3 动量预测，非正交循环？

OpenFOAM手册中对动量预测开关的建议是多相流情况时关闭，溃坝算例中的默认设置也默认是关闭的。从个人理解来看，动量预测可以提高求解的稳定性，但是会增加计算资源的消耗。

手册对于非正交循环次数究竟应该设置为多少则比较暧昧。按照手册的说法，非正交循环是为了校正在压力泊松方程中由于网格非正交引起的面法向梯度计算误差,但其校正的原理究竟是什么我暂时不清楚。根据对代码的分析和实际测试表明，计算中一个时间步内总的压力修正方程求解次数等于设置的(非正交循环次数+1)×设置的压力泊松方程求解次数(nCorrectors)，所谓非正交循环只是重复求解了一次压力泊松方程而已，这样的重复求解究竟如何实现非正交项的修正呢？另外，手册上建议在瞬态时将非正交修正次数设为1，稳态时将非正交修正次数设为0，但溃坝的瞬态算例却依然将非正交修正次数设置为0，这又是矛盾的。

根据我们的结果很难说明动量预测方程和非正交循环次数究竟应该设置为多少，或许各位读者可以使用其他比较好定量分析的算例对这个问题进行探究。根据我的经验，开启动量预测和增加非正交修正次数可以增强计算是稳定性（有时如果不开启非正交修正会导致计算误差累累积而导致错误结果），如果计算资源允许，可以开启这两个选项。

4 paraView后处理分析

本文在制作上面的动图时用到了paraview中的几个技巧：
绘制界面：
在二维的情况下，需要将counter滤镜的显示方式设置为Wireframe才可以与云图叠加显示。另外，还可以使用Threshold滤镜显示所有包含液体的区域：
添加标签：
添加时间标签：AnnotateTimeFilter滤镜
添加普通标签：Source菜单-Text