InterFoam

InterFoam Solver for 2 incompressible, isothermal immiscible fluids using a VOF

   (volume of fluid) phase-fraction based interface capturing approach,
   with optional mesh motion and mesh topology changes including adaptive
   re-meshing.

1 Related information in the web

1.1 Online discussion

• VOF method
• About interFoam solver
• How to set the two parameters of solver 'interFoam'?
• Inlet in interFoam (recovered via web.archive.org)
• Question about interface flow and adaptive mesh - In this thread Henry Weller explains how interFoam is different from the CICSAM method.

2 Equations

The solver solves the Navier Stokes equations for two incompressible, isothermal immiscible fluids. That means that the material properties are constant in the region filled by one of the two fluid except at the interphase.

2.1 Continuity Equation

The constant-density continuity equation is

$\frac{\partial {u}_j}{\partial x_j} = 0$

(1)

2.2 Momentum Equation

$\frac{ \partial (\rho {u}_i)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_j u_i \right) = - \frac{\partial p} {\partial{x_i}} + \frac{\partial}{\partial x_j} \left( \tau_{ij} + \tau_{t_{ij}} \right) + \rho g_i + \sigma_i,$

(2)

$u$ represent the velocity, $g_i$ the gravitational acceleration, $p$ the pressure and $\tau_{ij}$ and $\tau_{t_{ij}}$ are the viscose and turbulent stresses. $\sigma_i$ is the surface tension.

The density $\rho$ is defined as follows:

$\rho = \alpha \rho_1 + (1 - \alpha) \rho_2$

(3)

$\alpha$ is 1 inside fluid 1 with the density $\rho_1$ and 0 inside fluid 2 with the density $\rho_2$. At the interphase between the two fluids $\alpha$ varies between 0 and 1.

2.3 Equation for the interphase

In order to know where the interphase between the two fluids is, an additional equation for $\alpha$ has to be solved.

$\frac{\partial \alpha}{\partial t} + \frac{\partial ( \alpha {u}_j )}{\partial x_j} = 0$

(4)

The equation can be seen as the conservation of the mixture components along the path of a fluid parcel.

3 Source Code

\\OpenFOAM 6

3.1 interFoam.C

 
 
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "CMULES.H"
#include "EulerDdtScheme.H"
#include "localEulerDdtScheme.H"
#include "CrankNicolsonDdtScheme.H"
#include "subCycle.H"
#include "immiscibleIncompressibleTwoPhaseMixture.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "CorrectPhi.H"
#include "fvcSmooth.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
int main(int argc, char *argv[])
{
    #include "postProcess.H"
 
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
    #include "initContinuityErrs.H"
    #include "createDyMControls.H"
    #include "createFields.H"
    #include "createAlphaFluxes.H"
    #include "initCorrectPhi.H"
    #include "createUfIfPresent.H"
 
    turbulence->validate();
 
    if (!LTS)
    {
        #include "CourantNo.H"
        #include "setInitialDeltaT.H"
    }
 
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    Info<< "\nStarting time loop\n" << endl;
 
    while (runTime.run())
    {
        #include "readDyMControls.H"
 
        if (LTS)
        {
            #include "setRDeltaT.H"
        }
        else
        {
            #include "CourantNo.H"
            #include "alphaCourantNo.H"
            #include "setDeltaT.H"
        }
 
        runTime++;
 
        Info<< "Time = " << runTime.timeName() << nl << endl;
 
        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            if (pimple.firstIter() || moveMeshOuterCorrectors)
            {
                mesh.update();
 
                if (mesh.changing())
                {
                    // Do not apply previous time-step mesh compression flux
                    // if the mesh topology changed
                    if (mesh.topoChanging())
                    {
                        talphaPhi1Corr0.clear();
                    }
 
                    gh = (g & mesh.C()) - ghRef;
                    ghf = (g & mesh.Cf()) - ghRef;
 
                    MRF.update();
 
                    if (correctPhi)
                    {
                        // Calculate absolute flux
                        // from the mapped surface velocity
                        phi = mesh.Sf() & Uf();
 
                        #include "correctPhi.H"
 
                        // Make the flux relative to the mesh motion
                        fvc::makeRelative(phi, U);
 
                        mixture.correct();
                    }
 
                    if (checkMeshCourantNo)
                    {
                        #include "meshCourantNo.H"
                    }
                }
            }
 
            #include "alphaControls.H"
            #include "alphaEqnSubCycle.H"
 
            mixture.correct();
 
            #include "UEqn.H"
 
            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }
 
            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }
 
        runTime.write();
 
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
 
    Info<< "End\n" << endl;
 
    return 0;
}
 
 
// ************************************************************************* //