Difference between revisions of "InterFoam"

From OpenFOAMWiki
(→‎Online discussion: recovered link from web.archive.org)
Line 1: Line 1:
−
'''InterFoam''' is a solver for 2 incompressible fluids, which tracks the interface and includes the option of mesh motion.
+
'''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.
 +
 
  
 
==Related information in the web==
 
==Related information in the web==
Line 10: Line 14:
 
*[http://www.cfd-online.com/Forums/openfoam-solving/59643-question-about-interface-flow-adaptive-mesh.html Question about interface flow and adaptive mesh] - In this thread Henry Weller explains how interFoam is different from the CICSAM method.
 
*[http://www.cfd-online.com/Forums/openfoam-solving/59643-question-about-interface-flow-adaptive-mesh.html Question about interface flow and adaptive mesh] - In this thread Henry Weller explains how interFoam is different from the CICSAM method.
  
−
// Version Openfoam 1.5
+
==Equations==
−
<cpp>
+
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.
 +
 
 +
===Continuity Equation===
 +
The constant-density continuity equation is
 +
 
 +
<table width="70%"><tr><td>
 +
<center><math>
 +
 
 +
  \frac{\partial {u}_j}{\partial x_j} = 0
 +
 
 +
</math></center>
 +
</td><td width="5%">(1)</td></tr></table>
 +
 
 +
===Momentum Equation===
 +
 
 +
 
 +
<table width="70%"><tr><td>
 +
<center><math>
 +
 
 +
    \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,
 +
 
 +
</math></center>
 +
</td><td width="5%">(2)</td></tr></table>
 +
 
 +
<math> u </math> represent the velocity, <math> g_i </math> the gravitational acceleration, <math> p </math> the pressure and
 +
<math> \tau_{ij}  </math> and <math> \tau_{t_{ij}}  </math> are the viscose and turbulent stresses.
 +
 
 +
 
 +
The density <math> \rho </math> is defined as follows:
 +
 
 +
<table width="70%"><tr><td>
 +
<center><math>
 +
 
 +
  \rho = \alpha  \rho_1 + (1 - \alpha)  \rho_2
 +
 
 +
</math></center>
 +
</td><td width="5%">(3)</td></tr></table>
 +
 
 +
<math> \alpha </math> is 1 inside fluid 1 with the density  <math> \rho_1 </math> and 0 inside fluid 2 with the density <math> \rho_2  </math>. At the interphase between
 +
the two fluids <math> \alpha </math> varies between 0 and 1.
 +
 
 +
===Equation for the interphase ===
 +
 
 +
In order to know where the interphase between the two fluids is, an additional equation for <math> \alpha </math> has to be solved.
 +
 
 +
==Source Code==
 +
\\OpenFOAM 6
 +
 
 +
===interFoam.C===
 +
 
 +
 
 +
<br><cpp>
 +
 
 
#include "fvCFD.H"
 
#include "fvCFD.H"
−
#include "MULES.H"
+
#include "dynamicFvMesh.H"
 +
#include "CMULES.H"
 +
#include "EulerDdtScheme.H"
 +
#include "localEulerDdtScheme.H"
 +
#include "CrankNicolsonDdtScheme.H"
 
#include "subCycle.H"
 
#include "subCycle.H"
−
#include "interfaceProperties.H"
+
#include "immiscibleIncompressibleTwoPhaseMixture.H"
−
#include "twoPhaseMixture.H"
+
#include "turbulentTransportModel.H"
 +
#include "pimpleControl.H"
 +
#include "fvOptions.H"
 +
#include "CorrectPhi.H"
 +
#include "fvcSmooth.H"
  
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Line 22: Line 89:
 
int main(int argc, char *argv[])
 
int main(int argc, char *argv[])
 
{
 
{
−
// these header files contain source code for common tasks that are used in numerous applications.
+
     #include "postProcess.H"
−
     #include "setRootCase.H"// set the run directory
+
 
−
    #include "createTime.H" // rootPath, caseName
+
     #include "setRootCaseLists.H"
−
     #include "createMesh.H" // timeName
+
     #include "createTime.H"
−
     #include "readEnvironmentalProperties.H"
+
     #include "createDynamicFvMesh.H"
−
     #include "readPISOControls.H"
+
 
     #include "initContinuityErrs.H"
 
     #include "initContinuityErrs.H"
 +
    #include "createDyMControls.H"
 
     #include "createFields.H"
 
     #include "createFields.H"
−
// Reading nGammaCorr, nGammaSubCycles
+
    #include "createAlphaFluxes.H"
−
// Creates pd, gamma, U, rho, rhoPhi(Mass flux)
+
    #include "initCorrectPhi.H"
−
// Creates and initialises the face-flux field phi.
+
    #include "createUfIfPresent.H"
−
// Reading transportProperties
+
−
// Construct interface from gamma distribution
+
  
−
     #include "readTimeControls.H"
+
     turbulence->validate();
−
    #include "correctPhi.H"
+
−
    #include "CourantNo.H"
+
−
    #include "setInitialDeltaT.H"
+
−
+
  
−
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
    if (!LTS)
 +
    {
 +
        #include "CourantNo.H"
 +
        #include "setInitialDeltaT.H"
 +
    }
  
 +
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     Info<< "\nStarting time loop\n" << endl;
  
−
// use the runTime object to control time stepping
 
 
     while (runTime.run())
 
     while (runTime.run())
 
     {
 
     {
−
         #include "readPISOControls.H"
+
         #include "readDyMControls.H"
−
        // Read nCorrectors, nNonOrthogonalCorrectors, momentumPredictor,
+
 
−
         //  fluxGradp, transSonic, nOuterCorrectors
+
         if (LTS)
−
         #include "readTimeControls.H"
+
         {
−
         // Read adjustTimeStep, maxCo, maxDeltaT
+
            #include "setRDeltaT.H"
−
         #include "CourantNo.H"
+
         }
−
        // Get the Mean and max Courant Numbers
+
         else
−
        #include "setDeltaT.H"
+
        {
−
         // Get the deltaT
+
            #include "CourantNo.H"
 +
            #include "alphaCourantNo.H"
 +
            #include "setDeltaT.H"
 +
         }
  
 
         runTime++;
 
         runTime++;
Line 63: Line 131:
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
  
−
         // Calculates the laminar average kinematic viscosity
+
         // --- Pressure-velocity PIMPLE corrector loop
−
         twoPhaseProperties.correct();
+
         while (pimple.loop())
 +
        {
 +
            if (pimple.firstIter() || moveMeshOuterCorrectors)
 +
            {
 +
                mesh.update();
  
−
        #include "gammaEqnSubCycle.H"
+
                if (mesh.changing())
−
        // Read nGammaCorr, nGammaSubCycles
+
                {
 +
                    // Do not apply previous time-step mesh compression flux
 +
                    // if the mesh topology changed
 +
                    if (mesh.topoChanging())
 +
                    {
 +
                        talphaPhi1Corr0.clear();
 +
                    }
  
−
        #include "UEqn.H"
+
                    gh = (g & mesh.C()) - ghRef;
 +
                    ghf = (g & mesh.Cf()) - ghRef;
  
−
        // --- PISO loop
+
                    MRF.update();
−
        for (int corr=0; corr<nCorr; corr++)
+
−
        {
+
−
            #include "pEqn.H"
+
−
        }
+
  
−
        #include "continuityErrs.H"
+
                    if (correctPhi)
 +
                    {
 +
                        // Calculate absolute flux
 +
                        // from the mapped surface velocity
 +
                        phi = mesh.Sf() & Uf();
  
−
        p = pd + rho*gh;
+
                        #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();
 
         runTime.write();
Line 90: Line 201:
 
     Info<< "End\n" << endl;
 
     Info<< "End\n" << endl;
  
−
     return(0);
+
     return 0;
 
}
 
}
  
  
 
// ************************************************************************* //
 
// ************************************************************************* //
 +
</cpp><br>
  
−
</cpp>
 
  
 
[[Category:Multiphase flow solvers]]
 
[[Category:Multiphase flow solvers]]

Revision as of 10:57, 23 September 2018

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

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,
(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.


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.

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

Retrieved from "https://openfoamwiki.net/index.php?title=InterFoam&oldid=19320"
  • GNU Free Documentation License 1.3
  • Powered by MediaWiki
  • Powered by Semantic MediaWiki