数值模拟两相流求解器i n t e r F o a m的应用

张恩臻

（同济大学桥梁工程系，中国 上海 200092）

0 引言

计算流体动力学（Computational Fluid Dynamics，简称CFD）是通过计算机数值计算和图像显示，对包含有流体流动和热传导等相关物理现象的系统所做的分析。与实验研究相比，CFD计算具有速度快、成本低、资料完备等优点。

目前流行的CFD计算软件有Fluent、OpenFOAM等，其中penFOAM是Open Field Operation And Manipulation的英文缩写，很形象的概括了这款软件（严格的说，应该是开源的C++代码包）的特点：开源、直接对场操作、运算和处理。由于代码是开源的，所有代码均用户可见、可改、可重用，该软件具有非常大灵活性，几乎可以实现所有的用户自定义边界条件（时间相关、空间相关均可）、所有的来流条件，用户可以任意修改的N-S方程表达式（包括添加源项，添加新的场量，如温度），支持几乎所有的网格形式（包括结构、非结构网格）。

两相流的数值模拟是近年来的CFD领域研究的热点和前沿课题之一，目前流行的CFD计算软件都有各自的两相流求解器，interFoam是OpenFOAM中最基础最完善的两相流求解器。

1 interFoam及其特点

OpenFOAM有许多两相流求解器，interFoam是其中之一，用来求解不掺混的两相流 （掺混指的是一种流体分布在另一种流体当中，如气泡在液体中的运动）。InterFoam求解器中对于两相流自由界面的捕捉即是使用的流体体积法（VOF），而对于空间离散则是采用的有限体积法（FVM）。

interFoam是当下众多用于模拟两相流动是数值方法和代码中的一个流体体积法(VOF)数值工具，其正越来越受到更广泛的关注和应用，Trujillo.M.F[1]等(2011)的文章中也说明了这一点。这个VOF求解器最早是被Ubbink用在FOAM框架中的，其从最初的版本到现在经历了不断的改进。其如今的代码格式是OpenFOAM软件下的C++库中的一部分，通常用来对一系列与计算方法有关的偏微分方程作有限体积离散。C++语言面向对象的技术特点使其能与众多同类数学工具相竞争，也使得在其拓展和改进时，对高阶矩阵的操作更加容易。除此之外，还有其他一些优良的特点使得这一工具更具竞争力，如很容易实现并行、支持前处理和后处理工具、错误检查机制、可选的时间和空间离散化方法等等。

2 数值计算方法

N-S方程是公认的适用于所有流体的普适性偏微分方程，时间控制的三维可压缩牛顿流体运动及热量传递控制方程为：

方程存在极少的解析解，而通常工程问题应用中都是以取得偏微分方程的数值解为目标。大多数情况下计算N-S方程的精确解所需的计算量是巨大的，对工程问题是不可行的。这时候就需要一些简化的假设及模型使实际计算问题变得经济上可行。在interFoam求解其中假设流体是非粘性的，采用的控制方程为欧拉方程。

OpenFOAM采用有限体积法对偏微分方程进行空间离散。计算区域被分解成各个控制体，即离散的体积元。偏微分方程的通量形式在不同体积元之间相互联系。因为通常情况下流出体积元的通量是等于流入其的通量，所以一般有限体积法是守恒的。所涉及的量值（速度、压力等）是计算在体积元的质心的。相邻质心之间的数值可以通过差值来确定。

interFoam求解器是使用流体体积法（VOF）来进行界面捕捉的。其使用一个称为相参数(phase fraction,同上节的体积分数F)的特性数来捕捉自由流体界面。在OpenFOAM 1.5dev版本中，这个相参数是用字母γ来表示的，更新的版本则是用α来表示。这个参数是通过流场中对流来运输的，参数输运方程如下：

在只含有气体的体积元中此参数值为0，在充满液体的体积元中此参数值为1。在含有气体和液体接触的自由界面的体积元中，此参数值介于0和1之间。

interFoam中所使用的VOF方法有人为分散水汽接触界面的趋向，保持接触面形状清晰的一个方法是在垂直于接触面出施加人工压力梯度。在OpenFOAM中，人工体积力梯度可以通过cGamma设置来进行控制。

3 interFoam的应用

表1对interFoam求解器最新的应用实例做了一个列举。

研究人员 求解器应用实例Raach[7]等(2011)在interFoam框架中使用了一种能量方程。求解器在没有使用表面张力模型的情况下模拟了热传递的过程。将相速度和波峰高值与实验结果做了对比。Roisman[8]等(2010) 模拟了多孔界面的水滴效应，并于实验中所得的水花溅射结果做了大量对比。Saha和Mitra[9](2009)考虑了动态和静态接触角的同时，计算了微流体毛细流动 （microfluidic capillary flow），并与实验结果相比较。Saha[10]等(2009)使用了宾汉塑性流体模型。对其在渠道流的速度剖面计算所得解析解进行了进一步论证。证明了其在v形漏斗实验中得到的稳定流动解。Trujillo[12]等(2011)在interFoam中采用了一种能量方程，用以计算水滴效应下的热传递。证明了Soriano等所得实验结果中的花冠结构（crown formation），并支持了Trujillo实验中的温度结果。Deshpande[13]等(2012) 与他们一个下冲喷射流实验中所得速度剖面和液体碎片剖面作了对比。

Jonas Andersson(2011)曾利用OpenFOAM软件中的interFoam求解器对波浪及波浪引起的应力进行了数值模拟，并与实测数据进行了对比以验证此求解器的准确度。

blockMesh是OpenFOAM附带的工具之一，可以用来生成简单几何形状的网格。同OpenFOAM其他工具一样，blockMesh也是通过路径文件夹(dictionary)控制的，这使得参数化网格更加方便。输入若干关键参数便可以生成复杂网格，在原来参数上改动少许数据，则可以生成完全不同的复杂网格。这让blockMesh工具非常适合用于有几何相似性的一族网格的分析研究，比如参数最优化对网格的影响的研究。

ParaView是一个开源的数据分析及可视化后处理的工具。通过语句分析接口 （Message Parsing Interface，MPI）ParaView可以并行运行，使得大量数据可以快速分析。同其他大部分后处理软件一样，ParaView可以交互地3D化地运行。

Suraj S Deshpande[14]等(2012)曾对OpenFOAM 中的interFoam求解器在求解两相流的性能方面做过评估。通过不同的确定性指标测试，他们对此求解器做了一个全面的整体评估，这些测试内容包括：（1）（动力学）纯对流确认性测试；（2）高韦伯数限制的动力学性能；（3）表面压力驱动流动的动力学性能。相对应于（1），在此测试中，interFoam求解器的性能表现能够与已有的代数算法——流体体积法(VOF)相比对；然而，值得注意的是，其同那些几何重建算法(geometric reconstruction scheme)还是不可比的。 在（2）中，对于高密度比(θ~103)惯性驱动流动的模拟，其余理论值和实验结果表现出了高度的一致性。在（3）中，表面张力的影响十分重要，压力表面张力演变的连续性和曲率的准确性是也是十分重要的，这也与Francois等(2006)的结论相符。与此同时，他们也作了其他的一些测试，主要结论有：（a）interFoam的算法能够保证压力和表面张力变化的连续性；（b）此求解器计算的曲率值与理论值有微少差异(10%左右)，而且在此过程中略有变动。为了减少假性流动的扰动影响，他们在Galusinski和Vigneauxd(2008)的研究上更进一步，提出来一下interFoam毛细流动模拟稳定性的准则：

Δt≤max(10τμ,0.1τρ)（3）

其中 τμ=μΔx/σ，τρ=；他们还对于雾化有关的一些毛细流动做了模拟，得到的结果与已有文献中的数据也较为符合。

［1］Trujillo M F,Alvarado J,Gehring E,et al.Numerical simulations and experimental characterization of heat transfer from a periodic impingement of droplets[J].Journal of Heat Transfer,2011,133(12):122201.

［2］Berberovic＇E,van Hinsberg N P,Jakirlic＇S,et al.Drop impact onto a liquid layer of finite thickness:Dynamics of the cavity evolution[J].Physical Review E,2009,79(3):036306.

［3］Maiwald A,Schwarze R.Numerical analysis of flow-induced gas entrainment in roll coating[J].Applied Mathematical Modelling,2011,35(7):3516-3526.

［4］Srinivasan V,Salazar A J,Saito K.Modeling the disintegration of modulated liquid jets using volume-of-fluid (VOF)methodology[J].Applied Mathematical Modelling,2011,35(8):3710-3730.

［5］Liu X,García M H.Three-dimensional numerical model with free water surface and mesh deformation for local sediment scour[J].Journal of waterway,port,coastal,and ocean engineering,2008,134(4):203-217.

［6］Ishimoto J,Sato F,Sato G.Computational prediction of the effect of microcavitation on an atomization mechanism in a gasoline injector nozzle[J].Journal of Engineering for Gas Turbines and Power,2010,132(8):082801.

［7］Raach H,Somasundaram S,Mitrovic J.Optimisation of turbulence wire spacing in falling films performed with OpenFOAM[J].Desalination,2011,267(1):118-119.

［8］Roisman I V,Weickgenannt C M,Lembach A N,et al.Drop impact close to a pore:experimental and numerical investigations[C]//ILASS—Europe,23rd annual conference on liquid atomization and spray systems,Brno,Czech Republic.2010.

［9］Ashish Saha A,Mitra S K.Effect of dynamic contact angle in a volume of fluid(VOF)model for a microfluidic capillary flow[J].Journal of colloid and interface science,2009,339(2):461-480.

［10］Saha A A,Mitra S K,Tweedie M,et al.Experimental and numerical investigation of capillary flow in SU8 and PDMS microchannels with integrated pillars[J].Microfluidics and nanofluidics,2009,7(4):451-465.

［11］Ashish Saha A,Mitra S K.Effect of dynamic contact angle in a volume of fluid (VOF)model for a microfluidic capillary flow [J].Journal of colloid and interface science,2009,339(2):461-480.

［12］Trujillo M F,Alvarado J,Gehring E,et al.Numerical simulations and experimental characterization of heat transfer from a periodic impingement of droplets[J].Journal of Heat Transfer,2011,133(12):122201.

［13］Deshpande S S,Trujillo M F,Wu X,et al.Computational and experimental characterization of a liquid jet plunging into a quiescent pool at shallow inclination[J].International Journal of Heat and Fluid Flow,2012,34:1-14.

［14］Deshpande S S,Anumolu L,Trujillo M F.Evaluating the performance of the twophase flow solver interFoam[J].Computational Science&Discovery,2012,5(1):014016.