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.
Contents
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.
1.2 Useful links
An overview about interfoam: https://www.cfdyna.com/Home/OpenFOAM/of_Tut_Web/of_Suites_Description.pdf
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
![]() | (1) |
2.2 Momentum Equation
![]() | (2) |
represent the velocity,
the gravitational acceleration,
the pressure and
and
are the viscose and turbulent stresses.
is the
surface tension.
The density is defined as follows:
![]() | (3) |
is 1 inside fluid 1 with the density
and 0 inside fluid 2 with the density
. At the interphase between
the two fluids
varies between 0 and 1.
The surface tension 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}](/images/math/5/7/e/57e585bfda1b664547e369e0a9edec73.png)
is the surface tension constant and
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 )](/images/math/c/c/d/ccd1f12ea1e16ef06fd54368fd82bb16.png)
2.3 Equation for the interphase
In order to know where the interphase between the two fluids is, an additional equation for has to be solved.
![]() | (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; } // ************************************************************************* //
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() - ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh) ) * 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); adjustPhi(phiHbyA, U, p_rgh); fvc::makeAbsolute(phiHbyA, U); } 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); 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 += scAlpha*mag(mesh.delta() & fvc::interpolate(symm(fvc::grad(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; } } 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:
"(i) Tests of kinematics showed excellent mass conservation and acceptable advection errors. For the notched disc in rotational flow, the rate of convergence was found to be close to unity. However, for a deforming globule, it was slightly improved (∼1.5). This is comparable to the existing methodologies, which have a similar degree of maturity. However, it is noticeably worse when compared with the latest algorithms which include subgrid level geometric reconstruction. In addition, a notable disparity is found between the local truncation and the global error.
(ii) The solver was tested on inertia-dominated flows and was found to perform extremely well. Even with modest grid resolutions, accurate physics was captured for the two unsteady problems. Similar observations have been presented in the literature for flows in this high Weber number regime.
(iii) In relation to surface tension-dominated flows, we performed several verification tests and showed that the momentum formulation in interFoam is indeed designed to yield a discrete balance between pressure gradient and surface tension. The performance of the solver is comparable with the face centered balanced force approach of Francois et al . An interesting finding in this work is regarding the amplification of magnitude spurious currents with the factor 1/ρ (section 3.3.1). It is shown that this factor amplifies an imbalance that originates from residuals in the solution of the pressure equation and not from any algorithmic inconsistencies.
(iv) Accuracy of curvature computation was also evaluated through a verification test involving a stationary 2D droplet. The finding is that the average curvature and the average pressure jump across the interface do not converge to the analytical prediction. Such a convergence to a slightly erroneous value has already been reported in the literature [32] for VoF methods. These curvature errors, however, did not create any noticeable issues, and good agreement with experiments was obtained.
(v) Regarding spurious currents, following [66], we showed that for interFoam, the growth of spurious currents can be controlled by choosing a time step.
(vi) Finally, regarding computations that involved atomization, interFoam was able to capture the physics even with modest levels of grid resolution. This was exemplified in the test cases corresponding to a 2D standing capillary wave, Rayleigh breakup of a laminar jet and capillary retraction of a liquid jet."
5 References
- ↑ 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.
- ↑ 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.
- ↑ 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.
- ↑ 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.
- ↑ 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.