# 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.


## 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 + f_{\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. $f_{\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.

The surface tension $f_{\sigma i},$ is modelled as continuum surface force (CSF)[1]. It is calculated as follows (see also [2]):

$f_{\sigma i} = \sigma \kappa \frac{\partial \alpha }{\partial x_i}$

$\sigma$ is the surface tension constant and $\kappa$ the curvature. The curvature can be approximated as follows [3] and [4]:

$\kappa = - \frac{\partial n_i}{\partial x_i} = - \frac{\partial}{\partial x_i} \left ( \frac{\partial \alpha / \partial x_i}{|\partial \alpha / \partial x_i|}\right )$

### 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())
{

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;
}

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

### 3.2 UEqn.H



MRF.correctBoundaryVelocity(U);

fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(rho, U)
==
fvOptions(rho, U)
);

UEqn.relax();

fvOptions.constrain(UEqn);

if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
mixture.surfaceTensionForce()
) * mesh.magSf()
)
);

fvOptions.correct(U);
}


### 3.3 pEqn.H


{
if (correctPhi)
{
rAU.ref() = 1.0/UEqn.A();
}
else
{
rAU = 1.0/UEqn.A();
}

surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
);
MRF.makeRelative(phiHbyA);

if (p_rgh.needReference())
{
fvc::makeRelative(phiHbyA, U);
fvc::makeAbsolute(phiHbyA, U);
}

surfaceScalarField phig
(
(
mixture.surfaceTensionForce()
)*rAUf*mesh.magSf()
);

phiHbyA += phig;

// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);

p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - p_rghEqn.flux();

p_rgh.relax();

U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}

#include "continuityErrs.H"

// Correct Uf if the mesh is moving
fvc::correctUf(Uf, U, phi);

// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);

p == p_rgh + rho*gh;

if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}

if (!correctPhi)
{
rAU.clear();
}
}



### 3.4 alphaEqn.H



{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");

// Set the off-centering coefficient according to ddt scheme
scalar ocCoeff = 0;
{
tmp<fv::ddtScheme<scalar>> tddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("ddt(alpha)")
)
);
const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();

if
(
isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
|| isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
)
{
ocCoeff = 0;
}
else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
{
if (nAlphaSubCycles > 1)
{
FatalErrorInFunction
<< "Sub-cycling is not supported "
"with the CrankNicolson ddt scheme"
<< exit(FatalError);
}

if
(
alphaRestart
|| mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
)
{
ocCoeff =
refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
.ocCoeff();
}
}
else
{
FatalErrorInFunction
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
}

// Set the time blending factor, 1 for Euler
scalar cnCoeff = 1.0/(1.0 + ocCoeff);

// Standard face-flux compression coefficient
surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));

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

// Add the optional shear compression contribution
if (scAlpha > 0)
{
phic +=
}

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;
}
}

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)
{
#include "alphaSuSp.H"

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)
// - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
//           + fvc::div(phiCN), alpha1)
==
Su + fvm::Sp(Sp + divU, alpha1)
);

alpha1Eqn.solve();

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> talphaPhi1UD(alpha1Eqn.flux());
alphaPhi10 = talphaPhi1UD();

if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
{
Info<< "Applying the previous iteration compression flux" << endl;
MULES::correct
(
geometricOneField(),
alpha1,
alphaPhi10,
talphaPhi1Corr0.ref(),
oneField(),
zeroField()
);

alphaPhi10 += talphaPhi1Corr0();
}

// Cache the upwind-flux
talphaPhi1Corr0 = talphaPhi1UD;

alpha2 = 1.0 - alpha1;

mixture.correct();
}

for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
#include "alphaSuSp.H"

surfaceScalarField phir(phic*mixture.nHatf());

tmp<surfaceScalarField> talphaPhi1Un
(
fvc::flux
(
phiCN(),
cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);

if (MULESCorr)
{
tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
volScalarField alpha10("alpha10", alpha1);

MULES::correct
(
geometricOneField(),
alpha1,
talphaPhi1Un(),
talphaPhi1Corr.ref(),
Sp,
(-Sp*alpha1)(),
oneField(),
zeroField()
);

// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
alphaPhi10 += talphaPhi1Corr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
alphaPhi10 += 0.5*talphaPhi1Corr();
}
}
else
{
alphaPhi10 = talphaPhi1Un;

MULES::explicitSolve
(
geometricOneField(),
alpha1,
phiCN,
alphaPhi10,
Sp,
(Su + divU*min(alpha1(), scalar(1)))(),
oneField(),
zeroField()
);
}

alpha2 = 1.0 - alpha1;

mixture.correct();
}

if (alphaApplyPrevCorr && MULESCorr)
{
talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
}
else
{
talphaPhi1Corr0.clear();
}

#include "rhofs.H"

if
(
word(mesh.ddtScheme("ddt(rho,U)"))
== fv::EulerDdtScheme<vector>::typeName
|| word(mesh.ddtScheme("ddt(rho,U)"))
== fv::localEulerDdtScheme<vector>::typeName
)
{
rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
}
else
{
if (ocCoeff > 0)
{
// Calculate the end-of-time-step alpha flux
alphaPhi10 =
(alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
}

// Calculate the end-of-time-step mass flux
rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f;
}

Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
}


## 4 Validation

In this section the validation performed in the literature a summarized.

[5] Performed an accurate validation of interfoam. They found that:

## 5 References

1. Brackbill, J. U., Douglas B. Kothe, and Charles Zemach. "A continuum method for modeling surface tension." Journal of computational physics 100.2 (1992): 335-354.
2. Heyns, Johan A., and Oliver F. Oxtoby. "Modelling surface tension dominated multiphase flows using the VOF approach." 6th European Conference on Computational Fluid Dynamics. 2014.
3. Brackbill, J. U., Douglas B. Kothe, and Charles Zemach. "A continuum method for modeling surface tension." Journal of computational physics 100.2 (1992): 335-354.
4. Heyns, Johan A., and Oliver F. Oxtoby. "Modelling surface tension dominated multiphase flows using the VOF approach." 6th European Conference on Computational Fluid Dynamics. 2014.
5. Deshpande, S. S., Anumolu, L., & Trujillo, M. F. (2012). Evaluating the performance of the two-phase flow solver interFoam. Computational science & discovery, 5(1), 014016.