interFoam是基于 VOF 模型的不可压、绝热、不可溶、不掺混(掺混指的是一种流体分布在另-种流体当中，如气泡在液体中的运动)两相界面捕获求解器。InterFoam 求解器中对于两相流自由界面的捕捉即是使用的流体体积法(VOF),而对于空间离散则是采用的有限体积法(FVM)。

本文参考了苏军伟的博客，对于我这种初学菜鸟来讲帮助很大，感兴趣的小伙伴可以前往查看。

代码解析

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    interFoam

Description
    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
    (volume of fluid) phase-fraction based interface capturing approach.

    The momentum and other fluid properties are of the "mixture" and a single
    momentum equation is solved.

    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.

    For a two-fluid approach see twoPhaseEulerFoam.

\*---------------------------------------------------------------------------*/
//在OpenFOAM的求解器中，涉及到的构建时间、组建矩阵、有限体积离散、组建网格、量纲设置等有限体积库相关的头文件。建议在自定义的求解器中必备此项。
#include "fvCFD.H"
//MULES算法求解相分数
#include "CMULES.H"
//欧拉离散格式
#include "EulerDdtScheme.H"
//局部欧拉离散格式
#include "localEulerDdtScheme.H"
//Crank-Nicolson隐式差分格式
#include "CrankNicolsonDdtScheme.H"
//亚循环，减少运算时间（相方程求解是显示的）
#include "subCycle.H"
//不可混不可压缩两相流
#include "immiscibleIncompressibleTwoPhaseMixture.H"
//湍流模型
#include "turbulentTransportModel.H"
//定义PIMPLE循环，使用PIMPLE循环必备头文件
#include "pimpleControl.H"
//一个可供选择的物理框架，用做控制方程的源项或约束
#include "fvOptions.H"
//通量修正
#include "CorrectPhi.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//进入主程序
int main(int argc, char *argv[])
{
    #include "postProcess.H"
    //根据输入参数argc和argv设置算例根目录rootcase,必备头文件
    #include "setRootCase.H"
    //创建时间对象 runTime 用于控制程序运行（非稳态求解器都需要此头文件）
    #include "createTime.H"
    //创建网格对象，必备头文件
    #include "createMesh.H"
    //创建控制字典
    #include "createControl.H"
    #include "createTimeControls.H"
    #include "createRDeltaT.H"
    //初始化连续性误差
    #include "initContinuityErrs.H"
    //位于求解器根目录下，用于场对象的输入与输出
    #include "createFields.H"
    //创建源项，无需源项可删除
    #include "createFvOptions.H"
    //通量修正
    #include "correctPhi.H"
    //进行湍流验证
    turbulence->validate();

    if (!LTS)  //LTS：局部时间格式（欧拉形式）
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setInitialDeltaT.H"
    }

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//显示开始时间循环
    Info<< "\nStarting time loop\n" << endl;
//进入时间循环
    while (runTime.run())     //与 while (runTime.loop())功能相同， 但需添加runTime++;
    {
        #include "readTimeControls.H"

        if (LTS)
        {
            #include "setRDeltaT.H"
        }
        else
        {
            #include "CourantNo.H"  //计算并输出平均库朗特数以及最大库朗特数
            #include "alphaCourantNo.H" //计算并输出混合后的平均库朗特数以及最大库朗特数
            #include "setDeltaT.H"  //重置timestep以保持恒定的最大courant数;时间步长的减少是立即的，但增加是避免不稳定振荡。
        }

        runTime++;
//显示当前运行时间
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            //读取相分数控制字典
            #include "alphaControls.H"
            //亚循环，计算相分数，如非LTS则调用`alphaEqn.H`
            #include "alphaEqnSubCycle.H"

            mixture.correct(); //进行更新修正
            //求解动量方程
            #include "UEqn.H"

            // --- Pressure corrector loop（压力修正）
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())     //判断是否需进行湍流迭代
            {
            //重新计算湍动黏性，也就是对turbulence->divDevReff(U)需要的量进行更新。
            //比如k－e模型中下面函数包括求解k－e方程和计算有效黏性系数。
                turbulence->correct();
            }
        }
//输出计算结果到文件中
        runTime.write();
//显示当前执行时间和cpu时间
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
//显示程序结束
    Info<< "End\n" << endl;
//返回 0
    return 0;
}

// ************************************************************************* //

动量方程："UEqn.H"

动量方程代码含义可参见： OpenFOAM中动量方程 UEqn 浅析。

具体计算方式可参见： 多材料线弹性固体的应力分析

//MRF:基于旋转坐标系下的N-S方程修正边界速度,不需要可忽略
MRF.correctBoundaryVelocity(U);
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
//该项包含两部分：1.分子扩散项 2.雷诺应力的偏分量(dev)的散度。
//雷诺应力主分量（hyd）的散度归结到了压力项中，这是大多数雷诺时均和大涡模型实现的一贯做法。
//因此压力比层流模型中多了一项,那就是雷诺应力的主分量，通常被称为湍动压力（动压：p_rgh）                   
      + turbulence->divDevRhoReff(rho, U)    
     ==
        fvOptions(rho, U)     //源项或约束，不需要可忽略
    );
    //松弛迭代,加速收敛
    UEqn.relax();
//更正边界条件，维持通量守恒
    fvOptions.constrain(UEqn);
//如果进行速度预测（on），则求解完整的动量方程得到预测速度
//如果不进行速度预测(off），则预测速度直接取当前已知时间步的速度
//momentumPredictor：动量预测求解开关，对多相流以及低雷诺数一般设置为off；
    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct   //显示重构网格中心面（光滑了整个分布情况），转换为体中心通量
            (
                (
                    mixture.surfaceTensionForce()     //表面张力项
                  - ghf*fvc::snGrad(rho)              //体积力项（垂直梯度）
                  - fvc::snGrad(p_rgh)                // `p_rgh:动压`;压力梯度项
                ) * mesh.magSf()                      //网格面矢量绝对值      
            )
        );
//保证动量守恒
        fvOptions.correct(U);
    }

压力方程 "pEqn.H"

{
    // rAU:在速度的最后一个解中，从矩阵中提取对角项并存储倒数
    volScalarField rAU("rAU", 1.0/UEqn.A());
    //转换为表面标量场
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
    //定义HbyA
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
    //将上述体积矢量场转化为面心标量场（参考有限体积法）,并保证速度通量的全局守恒，以确保压力方程有解
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::flux(HbyA)
      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
    );
    MRF.makeRelative(phiHbyA);//将绝对通量转化为相对于面的通量(参考动网格技术)
    //对方程进行修正，保证速度边界条件守恒
    adjustPhi(phiHbyA, U, p_rgh);
    surfaceScalarField phig
    (
        (
            mixture.surfaceTensionForce()
          - ghf*fvc::snGrad(rho)
        )*rAUf*mesh.magSf()
    );
    //更新面通量场
    phiHbyA += phig;

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
    // 非正交压力修正循环 ，根据fvSolution字典文件中设定值n,求解以下方程n-1次 
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );
        //只有求解区域所有的压力边界都为第二类边界条件时，上面的值才会用到。
        //如果有第一类边界条件，压力参考值为该点处边界值。
        //对于不可压缩流动压力值为相对值，上面的参考值的大小对结果无影响，除非存在边界压力。
        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
		//通过查询system/fvSolution下的PIMPLE更新压力参考网格并重新设定参考值
        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
        //在最后一次非正交内循环中，使用最新压力校正通量
        if (pimple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - p_rghEqn.flux();
			//若外迭代次数为n,则压力场的松弛仅在n-1次外循环前进行
			//若外迭代次数为1，这里使用piso算法，压力不进行亚松驰
            p_rgh.relax();
            //校正速度，满足边界条件(主要针对第二类边界条件）
            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }
    //计算连续性方程误差
    #include "continuityErrs.H"
   // 计算总压
    p == p_rgh + rho*gh;
   //更新动压
    if (p_rgh.needReference())
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );
        p_rgh = p - rho*gh;
    }
}

相分数求解头文件：alphaEqn.H

{
   //制定两种离散格式
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");
   //制定时间离散格式
    tmp<fv::ddtScheme<scalar>> ddtAlpha
    (
        fv::ddtScheme<scalar>::New
        (
            mesh,
            mesh.ddtScheme("ddt(alpha)")
        )
    );
    //设置偏心系数
    // Set the off-centering coefficient according to ddt scheme
    scalar ocCoeff = 0;
    if
    (
        isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
     || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
    )
    {
        ocCoeff = 0;     //一阶离散该系数为零
    }
    else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))   //Crank-Nicolson 离散格式（二阶离散）
    {
        if (nAlphaSubCycles > 1)
        {   
            //子循环次数>1，报错
            FatalErrorInFunction
                << "Sub-cycling is not supported "
                   "with the CrankNicolson ddt scheme"
                << exit(FatalError);
        }
       //按照设定的时间离散格式求解该系数
        ocCoeff =
            refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()) //const:不可修改
           .ocCoeff();
    }
    else
    {
        FatalErrorInFunction
            << "Only Euler and CrankNicolson ddt schemes are supported"   //不支持其他输出格式
            << exit(FatalError);
    }
    //规定标量cnCoeff系数值
    scalar cnCoeff = 1.0/(1.0 + ocCoeff);

    // Standard face-flux compression coefficient
    surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));  //`mixture.cAlpha`:压缩因子

    // Add the optional isotropic（各向同性） compression contribution
    if (icAlpha > 0)
    {
        phic *= (1.0 - icAlpha);
        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
    }

    surfaceScalarField::Boundary& phicBf =
        phic.boundaryFieldRef();

    // Do not compress interface at non-coupled boundary faces（不压缩非耦合界面）
    // (inlets, outlets etc.)
    forAll(phic.boundaryField(), patchi)
    {
        fvsPatchScalarField& phicp = phicBf[patchi];

        if (!phicp.coupled())
        {
            phicp == 0;   //对于非耦合界面，该值为零
        }
    }
//Crank-Nicolson格式
    tmp<surfaceScalarField> phiCN(phi);

    // Calculate the Crank-Nicolson off-centred volumetric flux
    //体积流量
    if (ocCoeff > 0)
    {
        phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
    }
  
    if (MULESCorr) 
    {
        fvScalarMatrix alpha1Eqn
        (
            (
                LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            )
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
        );
        //求解相分数值
        alpha1Eqn.solve();
       //显示相1体积分数的极大值和极小值
        Info<< "Phase-1 volume fraction = "
            << alpha1.weightedAverage(mesh.Vsc()).value()
            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
            << endl;
//迎风格式
        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
        alphaPhi = talphaPhiUD();

        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
        {
            Info<< "Applying the previous iteration compression flux" << endl;
            //半隐式算法
            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);

            alphaPhi += talphaPhiCorr0();
        }

        // Cache the upwind-flux
        talphaPhiCorr0 = talphaPhiUD;

        alpha2 = 1.0 - alpha1;

        mixture.correct();
    }


    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
        surfaceScalarField phir(phic*mixture.nHatf()); //phir：
//总体通量
        tmp<surfaceScalarField> talphaPhiUn  
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
               -fvc::flux(-phir, alpha2, alpharScheme),   //方程6
                alpha1,
                alpharScheme
            )
        );

        // Calculate the Crank-Nicolson off-centred alpha flux
        //相分数流量
        if (ocCoeff > 0)
        {
            talphaPhiUn =
                cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
        }
//半隐式离散
        if (MULESCorr)
        {
            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
            volScalarField alpha10("alpha10", alpha1);
            MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0);

            // Under-relax the correction for all but the 1st corrector
            if (aCorr == 0)
            {
                alphaPhi += talphaPhiCorr(); 
            }
            else
            {
                alpha1 = 0.5*alpha1 + 0.5*alpha10;
                alphaPhi += 0.5*talphaPhiCorr();
            }
        }
        else
        {
            alphaPhi = talphaPhiUn;    //一阶迎风格式
            //全显性求解
            MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
        }

        alpha2 = 1.0 - alpha1;
//
        mixture.correct();
    }

    if (alphaApplyPrevCorr && MULESCorr)
    {
        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
    }

    if
    (
        word(mesh.ddtScheme("ddt(rho,U)"))
     == fv::EulerDdtScheme<vector>::typeName
    )
    {
      //给出欧拉离散格式下的动量方程中`rhoPhi`的值
        rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;  
    }
    else
    {
        if (ocCoeff > 0)
        {
            // Calculate the end-of-time-step alpha flux
            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
        }

        // Calculate the end-of-time-step mass flux
        rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
    }
//显示相1的相分数极值
    Info<< "Phase-1 volume fraction = "
        << alpha1.weightedAverage(mesh.Vsc()).value()
        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
        << endl;
}

说明：
fvOptions具体含义可参考fvOptions 浅析