solvers::incompressibleVoF: New solver module for incompressible two-phase flow with VoF

executed with foamRun for single region simulations of foamMultiRun for
multi-region simulations.  Replaces interFoam and all the corresponding
tutorials have been updated and moved to tutorials/modules/incompressibleVoF.

Both incompressibleVoF and compressibleVoF solver modules are derived from the
common two-phase VoF base-class solvers::VoFSolver which handles the
complexities of VoF interface-compression, boundedness and conservation with
2nd-order schemes in space and time using the semi-implicit MULES limiter and
solution proceedure.  This maximises code re-use, improves readability and
simplifies maintenance.

Class
    Foam::solvers::incompressibleVoF

Description
    Solver module for 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.

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

    Either mixture or two-phase transport modelling may be selected.  In the
    mixture approach a single laminar, RAS or LES model is selected to model the
    momentum stress.  In the Euler-Euler two-phase approach separate laminar,
    RAS or LES selected models are selected for each of the phases.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient and steady simulations.

    Optional fvModels and fvConstraints are provided to enhance the simulation
    in many ways including adding various sources, Lagrangian
    particles, surface film etc. and constraining or limiting the solution.

SourceFiles
    incompressibleVoF.C

See also
    Foam::solvers::VoFSolver
    Foam::solvers::compressibleVoF

applications/solvers/modules/Allwmake

@@ -9,6 +9,8 @@ wmake $targetType incompressibleFluid
 wmake $targetType isothermalFluid
 wmake $targetType fluid
 wmake $targetType multicomponentFluid
+wmake $targetType VoFSolver
+incompressibleVoF/Allwmake $targetType $*
 compressibleVoF/Allwmake $targetType $*
 multiphaseEuler/Allwmake $targetType $*
 wmake $targetType solid

applications/solvers/modules/VoFSolver/Make/files

@@ -0,0 +1,7 @@
+setRDeltaT.C
+moveMesh.C
+alphaPredictor.C
+momentumPredictor.C
+VoFSolver.C
+
+LIB = $(FOAM_LIBBIN)/libVoFSolver

applications/solvers/modules/VoFSolver/Make/options

@@ -0,0 +1,22 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
+    -I$(LIB_SRC)/physicalProperties/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/VoF \
+    -I$(LIB_SRC)/twoPhaseModels/interfaceCompression/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lfluidSolver \
+    -lphysicalProperties \
+    -linterfaceCompression \
+    -linterfaceProperties \
+    -ltwoPhaseMixture \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvModels \
+    -lfvConstraints \
+    -lsampling

applications/solvers/modules/VoFSolver/VoFSolver.C

@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "VoFSolver.H"
+#include "localEulerDdtScheme.H"
+#include "CorrectPhi.H"
+#include "geometricZeroField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(VoFSolver, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::VoFSolver::correctCoNum()
+{
+    fluidSolver::correctCoNum(phi);
+
+    const scalarField sumPhi
+    (
+        interface.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    const scalar meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+
+    Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+        << " max: " << alphaCoNum << endl;
+}
+
+
+void Foam::solvers::VoFSolver::continuityErrors()
+{
+    fluidSolver::continuityErrors(phi);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::VoFSolver::VoFSolver
+(
+    fvMesh& mesh,
+    autoPtr<twoPhaseMixture> mixturePtr
+)
+:
+    fluidSolver(mesh),
+
+    mixture_(mixturePtr),
+    mixture(mixture_()),
+
+    alpha1(mixture.alpha1()),
+    alpha2(mixture.alpha2()),
+
+    alphaRestart
+    (
+        typeIOobject<surfaceScalarField>
+        (
+            IOobject::groupName("alphaPhi", alpha1.group()),
+            runTime.name(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ).headerOk()
+    ),
+
+    divAlphaName("div(phi,alpha)"),
+
+    U
+    (
+        IOobject
+        (
+            "U",
+            runTime.name(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    ),
+
+    phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.name(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(U) & mesh.Sf()
+    ),
+
+    interface(mixture, alpha1, alpha2, U),
+
+    buoyancy(mesh),
+
+    p_rgh(buoyancy.p_rgh),
+
+    rho(mixture.rho()),
+
+    rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.name(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    ),
+
+    alphaPhi1
+    (
+        IOobject
+        (
+            IOobject::groupName("alphaPhi", alpha1.group()),
+            runTime.name(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        phi*fvc::interpolate(alpha1)
+    ),
+
+    MRF(mesh)
+{
+    // Read the controls
+    read();
+
+    mesh.schemes().setFluxRequired(p_rgh.name());
+    mesh.schemes().setFluxRequired(alpha1.name());
+
+    if (alphaRestart)
+    {
+        Info << "Restarting alpha" << endl;
+    }
+
+    if (mesh.dynamic())
+    {
+        Info<< "Constructing face momentum Uf" << endl;
+
+        Uf = new surfaceVectorField
+        (
+            IOobject
+            (
+                "Uf",
+                runTime.name(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            fvc::interpolate(U)
+        );
+    }
+
+    if (transient())
+    {
+        correctCoNum();
+    }
+    else if (LTS)
+    {
+        Info<< "Using LTS" << endl;
+
+        trDeltaT = tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fv::localEulerDdt::rDeltaTName,
+                    runTime.name(),
+                    mesh,
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh,
+                dimensionedScalar(dimless/dimTime, 1),
+                extrapolatedCalculatedFvPatchScalarField::typeName
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::VoFSolver::~VoFSolver()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::solvers::VoFSolver::maxDeltaT() const
+{
+    const scalar maxAlphaCo
+    (
+        runTime.controlDict().lookup<scalar>("maxAlphaCo")
+    );
+
+    const scalar deltaT = fluidSolver::maxDeltaT();
+
+    if (alphaCoNum > small)
+    {
+        return min
+        (
+            deltaT,
+            maxAlphaCo/(alphaCoNum + small)*runTime.deltaTValue()
+        );
+    }
+    else
+    {
+        return deltaT;
+    }
+}
+
+
+void Foam::solvers::VoFSolver::preSolve()
+{
+    // Read the controls
+    read();
+
+    if (transient())
+    {
+        correctCoNum();
+    }
+    else if (LTS)
+    {
+        setRDeltaT();
+    }
+
+    fvModels().preUpdateMesh();
+
+    // Store divU from the previous mesh so that it can be mapped
+    // and used in correctPhi to ensure the corrected phi has the
+    // same divergence
+    if
+    (
+        correctPhi
+     && divergent()
+     && mesh.topoChanged())
+    {
+        // Construct and register divU for mapping
+        divU = new volScalarField
+        (
+            "divU0",
+            fvc::div(fvc::absolute(phi, U))
+        );
+    }
+
+    // Update the mesh for topology change, mesh to mesh mapping
+    const bool topoChanged = mesh.update();
+
+    // Do not apply previous time-step mesh compression flux
+    // if the mesh topology changed
+    if (topoChanged)
+    {
+        talphaPhi1Corr0.clear();
+    }
+}
+
+
+void Foam::solvers::VoFSolver::prePredictor()
+{
+    fvModels().correct();
+    alphaPredictor();
+}
+
+
+void Foam::solvers::VoFSolver::postCorrector()
+{}
+
+
+void Foam::solvers::VoFSolver::postSolve()
+{
+    divU.clear();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/VoFSolver/VoFSolver.H

@@ -0,0 +1,282 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+Class
+    Foam::solvers::VoFSolver
+
+Description
+    Solver module base-class for for 2 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.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Either mixture or two-phase transport modelling may be selected.  In the
+    mixture approach a single laminar, RAS or LES model is selected to model the
+    momentum stress.  In the Euler-Euler two-phase approach separate laminar,
+    RAS or LES selected models are selected for each of the phases.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient and steady simulations.
+
+    Optional fvModels and fvConstraints are provided to enhance the simulation
+    in many ways including adding various sources, Lagrangian
+    particles, surface film etc. and constraining or limiting the solution.
+
+SourceFiles
+    VoFSolver.C
+
+See also
+    Foam::solvers::fluidSolver
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VoFSolver_H
+#define VoFSolver_H
+
+#include "fluidSolver.H"
+#include "twoPhaseMixture.H"
+#include "interfaceProperties.H"
+#include "buoyancy.H"
+#include "pressureReference.H"
+#include "IOMRFZoneList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class VoFSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class VoFSolver
+:
+    public fluidSolver
+{
+
+protected:
+
+    // Phase properties
+
+        //- The compressible two-phase mixture
+        autoPtr<twoPhaseMixture> mixture_;
+
+        twoPhaseMixture& mixture;
+
+        //- Reference to the phase1-fraction
+        volScalarField& alpha1;
+
+        //- Reference to the phase2-fraction
+        volScalarField& alpha2;
+
+        //- Switch indicating if this is a restart
+        bool alphaRestart;
+
+        //- Name of the alpha convection scheme
+        const word divAlphaName;
+
+        //- Phase-fraction flux Courant number
+        scalar alphaCoNum;
+
+
+    // Kinematic properties
+
+        //- Velocity field
+        volVectorField U;
+
+        //- Volumetric-flux field
+        surfaceScalarField phi;
+
+
+    // Interface properties
+
+        interfaceProperties interface;
+
+
+    // Thermophysical properties
+
+        //- Buoyancy force
+        solvers::buoyancy buoyancy;
+
+        //- Reference to the buoyant pressure for buoyant cases
+        //  otherwise to the pressure
+        volScalarField& p_rgh;
+
+        //- Reference to the mixture continuity density field
+        const volScalarField& rho;
+
+
+    // Kinematic properties
+
+        //- Mass flux field
+        surfaceScalarField rhoPhi;
+
+        // Phase-1 volumetric flux
+        surfaceScalarField alphaPhi1;
+
+
+    // Optional models
+
+        //- MRF zone list
+        IOMRFZoneList MRF;
+
+
+    // Cached temporary fields
+
+        tmp<volScalarField> rAU;
+
+        //- MULES Correction
+        tmp<surfaceScalarField> talphaPhi1Corr0;
+
+        //- Pointer to the surface momentum field
+        //  used to recreate the flux after mesh-change
+        autoPtr<surfaceVectorField> Uf;
+
+        //- Pointer to the momentum divergence field
+        //  used in correctPhi to ensure the corrected phi has the
+        //  same divergence
+        autoPtr<volScalarField> divU;
+
+        //- Optional LTS reciprocal time-step field
+        tmp<volScalarField> trDeltaT;
+
+        //- Cached momentum matrix
+        //  shared between the momentum predictor and pressure corrector
+        tmp<fvVectorMatrix> tUEqn;
+
+
+
+private:
+
+    // Private Member Functions
+
+        //- Set rDeltaT for LTS
+        virtual void setRDeltaT();
+
+        //- Correct the cached Courant numbers
+        void correctCoNum();
+
+        //- Solve for the phase-fractions
+        void alphaSolve(const dictionary& alphaControls);
+
+        //- Solve for the phase-fractions
+        void alphaPredictor();
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Calculate and print the continuity errors
+        void continuityErrors();
+
+        //- Return the pressure reference
+        virtual const Foam::pressureReference& pressureReference() const = 0;
+
+        //- Is the flow divergent?
+        //  i.e. compressible or include phase-fraction sources
+        virtual bool divergent() = 0;
+
+        //- Calculate the alpha equation sources
+        virtual void alphaSuSp
+        (
+            tmp<volScalarField::Internal>& Su,
+            tmp<volScalarField::Internal>& Sp
+        ) = 0;
+
+        //- Return the momentum equation stress term
+        virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U) = 0;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("VoFSolver");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        VoFSolver(fvMesh& mesh, autoPtr<twoPhaseMixture>);
+
+        //- Disallow default bitwise copy construction
+        VoFSolver(const VoFSolver&) = delete;
+
+
+    //- Destructor
+    virtual ~VoFSolver();
+
+
+    // Member Functions
+
+        //- Return the current maximum time-step for stable solution
+        virtual scalar maxDeltaT() const;
+
+        //- Called at the start of the time-step, before the PIMPLE loop
+        virtual void preSolve();
+
+        //- Called at the start of the PIMPLE loop to move the mesh
+        virtual bool moveMesh();
+
+        //- Called at the start of the PIMPLE loop
+        virtual void prePredictor() = 0;
+
+        //- Construct and optionally solve the momentum equation
+        virtual void momentumPredictor();
+
+        //- Construct and solve the energy equation,
+        //  convert to temperature
+        //  and update thermophysical and transport properties
+        virtual void thermophysicalPredictor() = 0;
+
+        //- Construct and solve the pressure equation in the PISO loop
+        virtual void pressureCorrector() = 0;
+
+        //- Correct the momentum transport modelling
+        virtual void postCorrector();
+
+        //- Called after the PIMPLE loop at the end of the time-step
+        virtual void postSolve();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const VoFSolver&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/VoFSolver/alphaPredictor.C

@@ -1,5 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "VoFSolver.H"
+#include "subCycle.H"
+#include "interfaceCompression.H"
+#include "CMULES.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "fvcFlux.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::VoFSolver::alphaSolve(const dictionary& alphaControls)
 {
-    #include "alphaScheme.H"
+    const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
+
+    const label nAlphaCorr(alphaControls.lookup<label>("nAlphaCorr"));
+
+    const bool MULESCorr
+    (
+        alphaControls.lookupOrDefault<Switch>("MULESCorr", false)
+    );
+
+    // Apply the compression correction from the previous iteration
+    // Improves efficiency for steady-simulations but can only be applied
+    // once the alpha field is reasonably steady, i.e. fully developed
+    const bool alphaApplyPrevCorr
+    (
+        alphaControls.lookupOrDefault<Switch>("alphaApplyPrevCorr", false)
+    );
+
+    const word alphaScheme(mesh.schemes().div(divAlphaName)[1].wordToken());
+
+    ITstream compressionScheme
+    (
+        compressionSchemes.found(alphaScheme)
+      ? mesh.schemes().div(divAlphaName)
+      : ITstream
+        (
+            divAlphaName,
+            tokenList
+            {
+                word("Gauss"),
+                word("interfaceCompression"),
+                alphaScheme,
+                alphaControls.lookup<scalar>("cAlpha")
+            }
+        )
+    );
+
 
     // Set the off-centering coefficient according to ddt scheme
     scalar ocCoeff = 0;
@@ -66,7 +137,22 @@
         );
     }
 
-    #include "alphaSuSp.H"
+    tmp<volScalarField> divU;
+
+    if (divergent())
+    {
+        divU =
+        (
+            mesh.moving()
+          ? fvc::div(phiCN + mesh.phi())
+          : fvc::div(phiCN)
+        );
+    }
+
+    tmp<volScalarField::Internal> Su;
+    tmp<volScalarField::Internal> Sp;
+
+    alphaSuSp(Su, Sp);
 
     if (MULESCorr)
     {
@@ -122,7 +208,7 @@
 
         alpha2 = 1.0 - alpha1;
 
-        mixture.correct();
+        interface.correct();
     }
 
     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
@@ -216,7 +302,7 @@
         alpha2 = 1.0 - alpha1;
 
         // Correct only the mixture interface for the interface compression flux
-        mixture.interfaceProperties::correct();
+        interface.correct();
     }
 
     if (alphaApplyPrevCorr && MULESCorr)
@@ -251,3 +337,62 @@
         << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
         << endl;
 }
+
+
+void Foam::solvers::VoFSolver::alphaPredictor()
+{
+    const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
+
+    const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
+
+
+    if (nAlphaSubCycles > 1)
+    {
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+        tmp<volScalarField> trSubDeltaT;
+
+        if (LTS)
+        {
+            trSubDeltaT =
+                fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+        }
+
+        // Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1
+        tmp<surfaceScalarField> talphaPhi1
+        (
+            surfaceScalarField::New
+            (
+                "alphaPhi1",
+                mesh,
+                dimensionedScalar(alphaPhi1.dimensions(), 0)
+            )
+        );
+
+        List<volScalarField*> alphaPtrs({&alpha1, &alpha2});
+
+        for
+        (
+            subCycle<volScalarField, subCycleFields> alphaSubCycle
+            (
+                alphaPtrs,
+                nAlphaSubCycles
+            );
+            !(++alphaSubCycle).end();
+        )
+        {
+            alphaSolve(alphaControls);
+            talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
+        }
+
+        alphaPhi1 = talphaPhi1();
+    }
+    else
+    {
+        alphaSolve(alphaControls);
+    }
+
+    mixture.correct();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/VoFSolver/alphaSuSp.H

@@ -2,10 +2,10 @@ tmp<volScalarField> divU;
 tmp<volScalarField::Internal> Su;
 tmp<volScalarField::Internal> Sp;
 
-if (fvModels.addsSupToField(alpha1.name()))
+if (divergent())
 {
     // Phase change alpha1 source
-    const fvScalarMatrix alphaSup(fvModels.source(alpha1));
+    const fvScalarMatrix alphaSup(fvModels().source(alpha1));
 
     Su = alphaSup.Su();
     Sp = alphaSup.Sp();

applications/solvers/modules/VoFSolver/momentumPredictor.C

@@ -23,18 +23,20 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "compressibleVoF.H"
+#include "VoFSolver.H"
+#include "fvmDiv.H"
+#include "fvcSnGrad.H"
+#include "fvcReconstruct.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::solvers::compressibleVoF::momentumPredictor()
+void Foam::solvers::VoFSolver::momentumPredictor()
 {
     tUEqn =
     (
         fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      - fvm::Sp(contErr1() + contErr2(), U)
       + MRF.DDt(rho, U)
-      + momentumTransport.divDevTau(U)
+      + divDevTau(U)
      ==
         fvModels().source(rho, U)
     );
@@ -53,7 +55,7 @@ void Foam::solvers::compressibleVoF::momentumPredictor()
             fvc::reconstruct
             (
                 (
-                    mixture.surfaceTensionForce()
+                    interface.surfaceTensionForce()
                   - buoyancy.ghf*fvc::snGrad(rho)
                   - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()
@@ -61,8 +63,6 @@ void Foam::solvers::compressibleVoF::momentumPredictor()
         );
 
         fvConstraints().constrain(U);
-
-        K = 0.5*magSqr(U);
     }
 }
 

applications/solvers/modules/VoFSolver/moveMesh.C

@@ -23,16 +23,22 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "compressibleVoF.H"
+#include "VoFSolver.H"
 #include "CorrectPhi.H"
+#include "geometricZeroField.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::solvers::compressibleVoF::moveMesh()
+bool Foam::solvers::VoFSolver::moveMesh()
 {
     if (pimple.firstIter() || pimple.moveMeshOuterCorrectors())
     {
-        if (correctPhi && !divU.valid())
+        if
+        (
+            correctPhi
+         && divergent()
+         && !divU.valid()
+        )
         {
             // Construct and register divU for correctPhi
             divU = new volScalarField
@@ -45,13 +51,13 @@ bool Foam::solvers::compressibleVoF::moveMesh()
         // Move the mesh
         mesh.move();
 
-        if (correctPhi || mesh.topoChanged())
+        if (mesh.changing())
         {
             buoyancy.moveMesh();
 
             MRF.update();
 
-            if (correctPhi)
+            if (correctPhi || mesh.topoChanged())
             {
                 // Calculate absolute flux
                 // from the mapped surface velocity
@@ -61,26 +67,42 @@ bool Foam::solvers::compressibleVoF::moveMesh()
 
                 if (correctPhi)
                 {
-                    CorrectPhi
-                    (
-                        phi,
-                        U,
-                        p_rgh,
-                        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-                        divU(),
-                        pressureReference,
-                        pimple
-                    );
+                    if (divU.valid())
+                    {
+                        CorrectPhi
+                        (
+                            phi,
+                            U,
+                            p_rgh,
+                            surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+                            divU(),
+                            pressureReference(),
+                            pimple
+                        );
+                    }
+                    else
+                    {
+                        CorrectPhi
+                        (
+                            phi,
+                            U,
+                            p_rgh,
+                            surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+                            geometricZeroField(),
+                            pressureReference(),
+                            pimple
+                        );
+                    }
                 }
 
                 // Make the fluxes relative to the mesh motion
                 fvc::makeRelative(phi, U);
             }
 
-            mixture.correct();
-
             meshCourantNo();
 
+            interface.correct();
+
             divU.clear();
 
             return true;

applications/solvers/modules/VoFSolver/setRDeltaT.C

@@ -23,12 +23,14 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "compressibleVoF.H"
+#include "VoFSolver.H"
 #include "fvcSmooth.H"
+#include "fvcSurfaceIntegrate.H"
+#include "fvcAverage.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::solvers::compressibleVoF::setRDeltaT()
+void Foam::solvers::VoFSolver::setRDeltaT()
 {
     volScalarField& rDeltaT = trDeltaT.ref();
 

applications/solvers/modules/compressibleVoF/Make/files

@@ -1,7 +1,4 @@
-setRDeltaT.C
-moveMesh.C
-alphaPredictor.C
-momentumPredictor.C
+alphaSuSp.C
 thermophysicalPredictor.C
 pressureCorrector.C
 compressibleVoF.C

applications/solvers/modules/compressibleVoF/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
     -IcompressibleInterPhaseTransportModel/lnInclude \
     -IcompressibleInterPhaseThermophysicalTransportModel/lnInclude \
@@ -9,7 +10,6 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
@@ -18,14 +18,14 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
-    -lfluidSolver \
+    -lVoFSolver \
     -lcompressibleInterPhaseTransportModel \
     -lcompressibleInterPhaseThermophysicalTransportModel \
     -lphysicalProperties \
     -lfluidThermophysicalModels \
     -linterfaceCompression \
     -lcompressibleInterfaceProperties \
-    -limmiscibleCompressibleTwoPhaseMixture \
+    -lcompressibleTwoPhaseMixture \
     -lfiniteVolume \
     -lmeshTools \
     -lfvModels \

applications/solvers/modules/compressibleVoF/alphaPredictor.C

@@ -32,79 +32,58 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::solvers::compressibleVoF::alphaPredictor()
+void Foam::solvers::compressibleVoF::alphaSuSp
+(
+    const surfaceScalarField& phiCN,
+    tmp<volScalarField>& divU,
+    tmp<volScalarField::Internal>& Su,
+    tmp<volScalarField::Internal>& Sp
+)
 {
-    #include "alphaControls.H"
-
-    const volScalarField& rho1 = mixture.thermo1().rho();
-    const volScalarField& rho2 = mixture.thermo2().rho();
-
-    tmp<surfaceScalarField> talphaPhi1(alphaPhi1);
-
-    if (nAlphaSubCycles > 1)
-    {
-        dimensionedScalar totalDeltaT = runTime.deltaT();
-
-        tmp<volScalarField> trSubDeltaT;
-
-        if (LTS)
-        {
-            trSubDeltaT =
-                fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
-        }
-
-        // Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1
-        tmp<surfaceScalarField> talphaPhi1
-        (
-            surfaceScalarField::New
-            (
-                "alphaPhi1",
-                mesh,
-                dimensionedScalar(alphaPhi1.dimensions(), 0)
-            )
-        );
-
-        // Sub-cycle on both alpha1 and alpha2
-        List<volScalarField*> alphaPtrs({&alpha1, &alpha2});
-
-        for
-        (
-            subCycle<volScalarField, subCycleFields> alphaSubCycle
-            (
-                alphaPtrs,
-                nAlphaSubCycles
-            );
-            !(++alphaSubCycle).end();
-        )
-        {
-            #include "alphaEqn.H"
-            talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
-        }
-
-        alphaPhi1 = talphaPhi1();
-    }
-    else
-    {
-        #include "alphaEqn.H"
-    }
-
-    mixture.correct();
-
-    alphaRhoPhi1 = fvc::interpolate(rho1)*alphaPhi1;
-    alphaRhoPhi2 = fvc::interpolate(rho2)*(phi - alphaPhi1);
-
-    rhoPhi = alphaRhoPhi1 + alphaRhoPhi2;
-
-    contErr1 =
+    Sp = volScalarField::Internal::New
     (
-        fvc::ddt(alpha1, rho1)()() + fvc::div(alphaRhoPhi1)()()
-      - (fvModels().source(alpha1, rho1)&rho1)()
+        "Sp",
+        mesh,
+        dimensionedScalar(dgdt.dimensions(), 0)
     );
 
-    contErr2 =
+    Su = volScalarField::Internal::New
     (
-        fvc::ddt(alpha2, rho2)()() + fvc::div(alphaRhoPhi2)()()
-      - (fvModels().source(alpha2, rho2)&rho2)()
+        "Su",
+        mesh,
+        dimensionedScalar(dgdt.dimensions(), 0)
+    );
+
+    if (fvModels().addsSupToField(alpha1.name()))
+    {
+        // Phase change alpha1 source
+        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
+
+        Su = alphaSup.Su();
+        Sp = alphaSup.Sp();
+    }
+
+    volScalarField::Internal& SpRef = Sp.ref();
+    volScalarField::Internal& SuRef = Su.ref();
+
+    forAll(dgdt, celli)
+    {
+        if (dgdt[celli] > 0.0)
+        {
+            SpRef[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+            SuRef[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        }
+        else if (dgdt[celli] < 0.0)
+        {
+            SpRef[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+        }
+    }
+
+    divU =
+    (
+        mesh.moving()
+      ? fvc::div(phiCN + mesh.phi())
+      : fvc::div(phiCN)
     );
 }
 

applications/solvers/modules/compressibleVoF/alphaSuSp.C

@@ -23,31 +23,55 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "immiscibleCompressibleTwoPhaseMixture.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::immiscibleCompressibleTwoPhaseMixture::
-immiscibleCompressibleTwoPhaseMixture(const volVectorField& U)
-:
-    compressibleTwoPhaseMixture(U.mesh()),
-    interfaceProperties(*this, alpha1(), alpha2(), U)
-{}
-
+#include "compressibleVoF.H"
+#include "fvcDiv.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::immiscibleCompressibleTwoPhaseMixture::correct()
+void Foam::solvers::compressibleVoF::alphaSuSp
+(
+    tmp<volScalarField::Internal>& Su,
+    tmp<volScalarField::Internal>& Sp
+)
 {
-    compressibleTwoPhaseMixture::correct();
-}
+    Sp = volScalarField::Internal::New
+    (
+        "Sp",
+        mesh,
+        dimensionedScalar(dgdt.dimensions(), 0)
+    );
 
+    Su = volScalarField::Internal::New
+    (
+        "Su",
+        mesh,
+        dimensionedScalar(dgdt.dimensions(), 0)
+    );
 
-bool Foam::immiscibleCompressibleTwoPhaseMixture::read()
-{
-    return
-        compressibleTwoPhaseMixture::read()
-     && interfaceProperties::read();
+    if (fvModels().addsSupToField(alpha1.name()))
+    {
+        // Phase change alpha1 source
+        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
+
+        Su = alphaSup.Su();
+        Sp = alphaSup.Sp();
+    }
+
+    volScalarField::Internal& SpRef = Sp.ref();
+    volScalarField::Internal& SuRef = Su.ref();
+
+    forAll(dgdt, celli)
+    {
+        if (dgdt[celli] > 0.0)
+        {
+            SpRef[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+            SuRef[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        }
+        else if (dgdt[celli] < 0.0)
+        {
+            SpRef[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+        }
+    }
 }
 
 

applications/solvers/modules/compressibleVoF/alphaSuSp.H

@@ -1,51 +0,0 @@
-tmp<volScalarField::Internal> Sp
-(
-    volScalarField::Internal::New
-    (
-        "Sp",
-        mesh,
-        dimensionedScalar(dgdt.dimensions(), 0)
-    )
-);
-
-tmp<volScalarField::Internal> Su
-(
-    volScalarField::Internal::New
-    (
-        "Su",
-        mesh,
-        dimensionedScalar(dgdt.dimensions(), 0)
-    )
-);
-
-if (fvModels().addsSupToField(alpha1.name()))
-{
-    // Phase change alpha1 source
-    const fvScalarMatrix alphaSup(fvModels().source(alpha1));
-
-    Su = alphaSup.Su();
-    Sp = alphaSup.Sp();
-}
-
-volScalarField::Internal& SpRef = Sp.ref();
-volScalarField::Internal& SuRef = Su.ref();
-
-forAll(dgdt, celli)
-{
-    if (dgdt[celli] > 0.0)
-    {
-        SpRef[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-        SuRef[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-    }
-    else if (dgdt[celli] < 0.0)
-    {
-        SpRef[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
-    }
-}
-
-tmp<volScalarField> divU
-(
-    mesh.moving()
-  ? fvc::div(phiCN() + mesh.phi())
-  : fvc::div(phiCN())
-);

applications/solvers/modules/compressibleVoF/compressibleInterPhaseThermophysicalTransportModel/Make/options

@@ -5,7 +5,6 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
@@ -17,7 +16,7 @@ EXE_INC = \
 LIB_LIBS = \
     -lphysicalProperties \
     -lfluidThermophysicalModels \
-    -limmiscibleCompressibleTwoPhaseMixture \
+    -lcompressibleTwoPhaseMixture \
     -lmomentumTransportModels \
     -lcompressibleMomentumTransportModels \
     -lphaseCompressibleMomentumTransportModels \

applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/Make/options

@@ -4,7 +4,6 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
@@ -15,7 +14,7 @@ EXE_INC = \
 LIB_LIBS = \
     -lphysicalProperties \
     -lfluidThermophysicalModels \
-    -limmiscibleCompressibleTwoPhaseMixture \
+    -lcompressibleTwoPhaseMixture \
     -lmomentumTransportModels \
     -lcompressibleMomentumTransportModels \
     -lphaseCompressibleMomentumTransportModels \

applications/solvers/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H

@@ -45,7 +45,7 @@ SourceFiles
 #define compressibleInterPhaseTransportModel_H
 
 #include "twoPhaseMixture.H"
-#include "immiscibleCompressibleTwoPhaseMixture.H"
+#include "compressibleTwoPhaseMixture.H"
 #include "compressibleMomentumTransportModels.H"
 #include "phaseCompressibleMomentumTransportModel.H"
 

applications/solvers/modules/compressibleVoF/compressibleVoF.C

@@ -39,82 +39,19 @@ namespace solvers
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::solvers::compressibleVoF::correctCoNum()
-{
-    fluidSolver::correctCoNum(phi);
-
-    const scalarField sumPhi
-    (
-        mixture.nearInterface()().primitiveField()
-       *fvc::surfaceSum(mag(phi))().primitiveField()
-    );
-
-    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
-
-    const scalar meanAlphaCoNum =
-        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
-
-    Info<< "Interface Courant Number mean: " << meanAlphaCoNum
-        << " max: " << alphaCoNum << endl;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
 :
-    fluidSolver(mesh),
-
-    U
+    VoFSolver
     (
-        IOobject
-        (
-            "U",
-            runTime.name(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
+        mesh,
+        autoPtr<twoPhaseMixture>(new compressibleTwoPhaseMixture(mesh))
     ),
 
-    phi
-    (
-        IOobject
-        (
-            "phi",
-            runTime.name(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        linearInterpolate(U) & mesh.Sf()
-    ),
+    mixture(refCast<compressibleTwoPhaseMixture>(VoFSolver::mixture)),
 
-    mixture(U),
-
-    alpha1(mixture.alpha1()),
-    alpha2(mixture.alpha2()),
-
-    alphaRestart
-    (
-        typeIOobject<surfaceScalarField>
-        (
-            IOobject::groupName("alphaPhi", alpha1.group()),
-            runTime.name(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ).headerOk()
-    ),
-
-    buoyancy(mesh),
-
-    p_rgh(buoyancy.p_rgh),
-
-    rho(mixture.rho()),
+    p(mixture.p()),
 
     dgdt
     (
@@ -129,9 +66,9 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
         alpha1*fvc::div(phi)
     ),
 
-    pressureReference
+    pressureReference_
     (
-        mixture.p(),
+        p,
         p_rgh,
         pimple.dict(),
         false
@@ -144,32 +81,6 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
         mixture
     ),
 
-    rhoPhi
-    (
-        IOobject
-        (
-            "rhoPhi",
-            runTime.name(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        fvc::interpolate(rho)*phi
-    ),
-
-    alphaPhi1
-    (
-        IOobject
-        (
-            IOobject::groupName("alphaPhi", alpha1.group()),
-            runTime.name(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        phi*fvc::interpolate(alpha1)
-    ),
-
     alphaRhoPhi1
     (
         IOobject::groupName("alphaRhoPhi", alpha1.group()),
@@ -196,66 +107,11 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
         mixture
     ),
 
-    thermophysicalTransport(momentumTransport),
-
-    MRF(mesh)
+    thermophysicalTransport(momentumTransport)
 {
     // Read the controls
     read();
 
-    mesh.schemes().setFluxRequired(p_rgh.name());
-    mesh.schemes().setFluxRequired(alpha1.name());
-
-    if (alphaRestart)
-    {
-        Info << "Restarting alpha" << endl;
-    }
-
-    if (mesh.dynamic())
-    {
-        Info<< "Constructing face momentum Uf" << endl;
-
-        Uf = new surfaceVectorField
-        (
-            IOobject
-            (
-                "Uf",
-                runTime.name(),
-                mesh,
-                IOobject::READ_IF_PRESENT,
-                IOobject::AUTO_WRITE
-            ),
-            fvc::interpolate(U)
-        );
-    }
-
-    if (transient())
-    {
-        correctCoNum();
-    }
-    else if (LTS)
-    {
-        Info<< "Using LTS" << endl;
-
-        trDeltaT = tmp<volScalarField>
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    fv::localEulerDdt::rDeltaTName,
-                    runTime.name(),
-                    mesh,
-                    IOobject::READ_IF_PRESENT,
-                    IOobject::AUTO_WRITE
-                ),
-                mesh,
-                dimensionedScalar(dimless/dimTime, 1),
-                extrapolatedCalculatedFvPatchScalarField::typeName
-            )
-        );
-    }
-
     if (correctPhi)
     {
         rAU = new volScalarField
@@ -283,77 +139,53 @@ Foam::solvers::compressibleVoF::~compressibleVoF()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::scalar Foam::solvers::compressibleVoF::maxDeltaT() const
-{
-    const scalar maxAlphaCo
-    (
-        runTime.controlDict().lookup<scalar>("maxAlphaCo")
-    );
-
-    const scalar deltaT = fluidSolver::maxDeltaT();
-
-    if (alphaCoNum > small)
-    {
-        return min
-        (
-            deltaT,
-            maxAlphaCo/(alphaCoNum + small)*runTime.deltaTValue()
-        );
-    }
-    else
-    {
-        return deltaT;
-    }
-}
-
-
-void Foam::solvers::compressibleVoF::preSolve()
-{
-    // Read the controls
-    read();
-
-    if (transient())
-    {
-        correctCoNum();
-    }
-    else if (LTS)
-    {
-        setRDeltaT();
-    }
-
-    // Store divU from the previous mesh so that it can be mapped
-    // and used in correctPhi to ensure the corrected phi has the
-    // same divergence
-    if (correctPhi && mesh.topoChanged())
-    {
-        // Construct and register divU for mapping
-        divU = new volScalarField
-        (
-            "divU0",
-            fvc::div(fvc::absolute(phi, U))
-        );
-    }
-
-    fvModels().preUpdateMesh();
-
-    // Update the mesh for topology change, mesh to mesh mapping
-    mesh.update();
-}
-
-
 void Foam::solvers::compressibleVoF::prePredictor()
 {
-    fvModels().correct();
-    alphaPredictor();
+    VoFSolver::prePredictor();
+
+    const volScalarField& rho1 = mixture.thermo1().rho();
+    const volScalarField& rho2 = mixture.thermo2().rho();
+
+    alphaRhoPhi1 = fvc::interpolate(rho1)*alphaPhi1;
+    alphaRhoPhi2 = fvc::interpolate(rho2)*(phi - alphaPhi1);
+
+    rhoPhi = alphaRhoPhi1 + alphaRhoPhi2;
+
+    contErr1 =
+    (
+        fvc::ddt(alpha1, rho1)()() + fvc::div(alphaRhoPhi1)()()
+      - (fvModels().source(alpha1, rho1)&rho1)()
+    );
+
+    contErr2 =
+    (
+        fvc::ddt(alpha2, rho2)()() + fvc::div(alphaRhoPhi2)()()
+      - (fvModels().source(alpha2, rho2)&rho2)()
+    );
 
     if (pimple.predictTransport())
     {
         momentumTransport.predict();
+    }
+
+    if (pimple.predictTransport())
+    {
         thermophysicalTransport.predict();
     }
 }
 
 
+void Foam::solvers::compressibleVoF::momentumPredictor()
+{
+    VoFSolver::momentumPredictor();
+
+    if (pimple.momentumPredictor())
+    {
+        K = 0.5*magSqr(U);
+    }
+}
+
+
 void Foam::solvers::compressibleVoF::postCorrector()
 {
     if (pimple.correctTransport())
@@ -364,10 +196,4 @@ void Foam::solvers::compressibleVoF::postCorrector()
 }
 
 
-void Foam::solvers::compressibleVoF::postSolve()
-{
-    divU.clear();
-}
-
-
 // ************************************************************************* //

applications/solvers/modules/compressibleVoF/compressibleVoF.H

@@ -49,15 +49,16 @@ SourceFiles
     compressibleVoF.C
 
 See also
-    Foam::solvers::fluidSolver
-    Foam::solvers::incompressibleFluid
+    Foam::solvers::VoFSolver
+    Foam::solvers::incompressibleVoF
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef compressibleVoF_H
 #define compressibleVoF_H
 
-#include "fluidSolver.H"
+#include "VoFSolver.H"
+#include "compressibleTwoPhaseMixture.H"
 #include "compressibleInterPhaseTransportModel.H"
 #include "compressibleInterPhaseThermophysicalTransportModel.H"
 #include "buoyancy.H"
@@ -76,48 +77,21 @@ namespace solvers
 
 class compressibleVoF
 :
-    public fluidSolver
+    public VoFSolver
 {
 
 protected:
 
-    // Kinematic properties
-
-        //- Velocity field
-        volVectorField U;
-
-        //- Volumetric-flux field
-        surfaceScalarField phi;
-
-
     // Phase properties
 
         //- The compressible two-phase mixture
-        immiscibleCompressibleTwoPhaseMixture mixture;
-
-        //- Reference to the phase1-fraction
-        volScalarField& alpha1;
-
-        //- Reference to the phase2-fraction
-        volScalarField& alpha2;
-
-        //- Switch indicating if this is a restart
-        bool alphaRestart;
-
-        scalar alphaCoNum;
+        compressibleTwoPhaseMixture& mixture;
 
 
     // Thermophysical properties
 
-        //- Buoyancy force
-        solvers::buoyancy buoyancy;
-
-        //- Reference to the buoyant pressure for buoyant cases
-        //  otherwise to the pressure
-        volScalarField& p_rgh;
-
-        //- Reference to the mixture continuity density field
-        const volScalarField& rho;
+        //- Reference to the mixture static pressure field
+        volScalarField& p;
 
         //- Compressibility source
         volScalarField::Internal dgdt;
@@ -126,7 +100,7 @@ protected:
     // Pressure reference
 
         //- Pressure reference
-        Foam::pressureReference pressureReference;
+        Foam::pressureReference pressureReference_;
 
         //- Minimum pressure
         dimensionedScalar pMin;
@@ -134,12 +108,6 @@ protected:
 
     // Kinematic properties
 
-        //- Mass flux field
-        surfaceScalarField rhoPhi;
-
-        // Phase-1 volumetric flux
-        surfaceScalarField alphaPhi1;
-
         // Phase-1 mass-flux
         surfaceScalarField alphaRhoPhi1;
 
@@ -159,56 +127,50 @@ protected:
 
     // Thermophysical transport
 
+        //- Thermophysical transport model
         compressibleInterPhaseThermophysicalTransportModel
             thermophysicalTransport;
 
 
-    // Optional models
-
-        //- MRF zone list
-        IOMRFZoneList MRF;
-
-
     // Cached temporary fields
 
-        tmp<volScalarField> rAU;
-
+        //- Phase-1 continuity error
         tmp<volScalarField::Internal> contErr1;
 
+        //- Phase-2 continuity error
         tmp<volScalarField::Internal> contErr2;
 
-        //- MULES Correction
-        tmp<surfaceScalarField> talphaPhi1Corr0;
 
-        //- Pointer to the surface momentum field
-        //  used to recreate the flux after mesh-change
-        autoPtr<surfaceVectorField> Uf;
+protected:
 
-        //- Pointer to the momentum divergence field
-        //  used in correctPhi to ensure the corrected phi has the
-        //  same divergence
-        autoPtr<volScalarField> divU;
+    // Protected Member Functions
 
-        //- Optional LTS reciprocal time-step field
-        tmp<volScalarField> trDeltaT;
+        //- Return the pressure reference
+        virtual const Foam::pressureReference& pressureReference() const
+        {
+            return pressureReference_;
+        }
 
-        //- Cached momentum matrix
-        //  shared between the momentum predictor and pressure corrector
-        tmp<fvVectorMatrix> tUEqn;
+        //- Compressible flow is divergent
+        virtual bool divergent()
+        {
+            return true;
+        }
 
+        //- Calculate the alpha equation sources
+        virtual void alphaSuSp
+        (
+            tmp<volScalarField::Internal>& Su,
+            tmp<volScalarField::Internal>& Sp
+        );
 
-private:
-
-    // Private Member Functions
-
-        //- Set rDeltaT for LTS
-        virtual void setRDeltaT();
-
-        //- Correct the cached Courant numbers
-        void correctCoNum();
-
-        //- Solve for the phase-fractions
-        void alphaPredictor();
+        //- Return the momentum equation stress term
+        virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U)
+        {
+            return
+                momentumTransport.divDevTau(U)
+              - fvm::Sp(contErr1() + contErr2(), U);
+        }
 
 
 public:
@@ -232,15 +194,6 @@ public:
 
     // Member Functions
 
-        //- Return the current maximum time-step for stable solution
-        virtual scalar maxDeltaT() const;
-
-        //- Called at the start of the time-step, before the PIMPLE loop
-        virtual void preSolve();
-
-        //- Called at the start of the PIMPLE loop to move the mesh
-        virtual bool moveMesh();
-
         //- Called at the start of the PIMPLE loop
         virtual void prePredictor();
 
@@ -258,9 +211,6 @@ public:
         //- Correct the momentum and thermophysical transport modelling
         virtual void postCorrector();
 
-        //- Called after the PIMPLE loop at the end of the time-step
-        virtual void postSolve();
-
 
     // Member Operators
 

applications/solvers/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/Make/options

@@ -4,7 +4,6 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
@@ -14,7 +13,7 @@ EXE_INC = \
 LIB_LIBS = \
     -lphysicalProperties \
     -lfluidThermophysicalModels \
-    -limmiscibleCompressibleTwoPhaseMixture \
+    -lcompressibleTwoPhaseMixture \
     -lmomentumTransportModels \
     -lcompressibleMomentumTransportModels \
     -lphaseCompressibleMomentumTransportModels \

applications/solvers/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C

@@ -24,8 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "VoFTurbulenceDamping.H"
-#include "immiscibleCompressibleTwoPhaseMixture.H"
-#include "interfaceProperties.H"
+#include "compressibleTwoPhaseMixture.H"
 #include "compressibleMomentumTransportModel.H"
 #include "fvMatrix.H"
 #include "addToRunTimeSelectionTable.H"
@@ -66,12 +65,11 @@ Foam::fv::compressible::VoFTurbulenceDamping::VoFTurbulenceDamping
     delta_("delta", dimLength, dict),
     mixture_
     (
-        mesh.lookupObject<immiscibleCompressibleTwoPhaseMixture>
+        mesh.lookupObject<compressibleTwoPhaseMixture>
         (
             "phaseProperties"
         )
     ),
-    interface_(mixture_),
     momentumTransport_
     (
         mesh.lookupType<compressibleMomentumTransportModel>(phaseName_)
@@ -146,13 +144,13 @@ void Foam::fv::compressible::VoFTurbulenceDamping::addSup
 
     if (fieldName == "epsilon")
     {
-        eqn += interface_.fraction()*C2_*aRhoSqrnu*momentumTransport_.k()()
-           /pow4(delta_);
+        eqn += mixture_.interfaceFraction()
+           *C2_*aRhoSqrnu*momentumTransport_.k()()/pow4(delta_);
     }
     else if (fieldName == "omega")
     {
-        eqn += interface_.fraction()*beta_*aRhoSqrnu
-           /(sqr(betaStar_)*pow4(delta_));
+        eqn += mixture_.interfaceFraction()
+           *beta_*aRhoSqrnu/(sqr(betaStar_)*pow4(delta_));
     }
     else
     {
@@ -197,13 +195,13 @@ void Foam::fv::compressible::VoFTurbulenceDamping::addSup
 
     if (fieldName == IOobject::groupName("epsilon", phaseName_))
     {
-        eqn += interface_.fraction()*C2_*taRhoSqrnu*momentumTransport_.k()()
-           /pow4(delta_);
+        eqn += mixture_.interfaceFraction()
+           *C2_*taRhoSqrnu*momentumTransport_.k()()/pow4(delta_);
     }
     else if (fieldName == IOobject::groupName("omega", phaseName_))
     {
-        eqn += interface_.fraction()*beta_*taRhoSqrnu
-           /(sqr(betaStar_)*pow4(delta_));
+        eqn += mixture_.interfaceFraction()
+           *beta_*taRhoSqrnu/(sqr(betaStar_)*pow4(delta_));
     }
     else
     {

applications/solvers/modules/compressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H

@@ -81,8 +81,7 @@ SourceFiles
 namespace Foam
 {
 
-class immiscibleCompressibleTwoPhaseMixture;
-class interfaceProperties;
+class compressibleTwoPhaseMixture;
 class compressibleMomentumTransportModel;
 
 namespace fv
@@ -111,10 +110,7 @@ class VoFTurbulenceDamping
         dimensionedScalar delta_;
 
         //- Reference to the mixture properties
-        const immiscibleCompressibleTwoPhaseMixture& mixture_;
-
-        //- Reference to the interface properties
-        const interfaceProperties& interface_;
+        const compressibleTwoPhaseMixture& mixture_;
 
         //- Reference to the mixture momentumTransport model
         const compressibleMomentumTransportModel& momentumTransport_;

applications/solvers/modules/compressibleVoF/pressureCorrector.C

@@ -32,8 +32,6 @@ License
 
 void Foam::solvers::compressibleVoF::pressureCorrector()
 {
-    volScalarField& p = mixture.p();
-
     const volScalarField& rho1 = mixture.rho1();
     const volScalarField& rho2 = mixture.rho2();
 
@@ -70,7 +68,7 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
         surfaceScalarField phig
         (
             (
-                mixture.surfaceTensionForce()
+                interface.surfaceTensionForce()
               - buoyancy.ghf*fvc::snGrad(rho)
             )*rAUf*mesh.magSf()
         );

applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C

@@ -39,7 +39,6 @@ void Foam::solvers::compressibleVoF::thermophysicalPredictor()
     const fvScalarMatrix e2Source(fvModels().source(alpha2, rho2, e2));
 
     volScalarField& T = mixture.T();
-    const volScalarField& p = mixture.p();
 
     fvScalarMatrix TEqn
     (

applications/solvers/modules/incompressibleVoF/Allwclean

@@ -7,6 +7,5 @@ wclean libso fvModels/VoFTurbulenceDamping
 wclean libso fvModels/VoFCavitation
 
 wclean
-wclean interMixingFoam
 
 #------------------------------------------------------------------------------

applications/solvers/modules/incompressibleVoF/Allwmake

@@ -10,6 +10,5 @@ wmake $targetType fvModels/VoFTurbulenceDamping
 wmake $targetType fvModels/VoFCavitation
 
 wmake $targetType
-wmake $targetType interMixingFoam
 
 #------------------------------------------------------------------------------

applications/solvers/modules/incompressibleVoF/Make/files

@@ -0,0 +1,5 @@
+alphaSuSp.C
+pressureCorrector.C
+incompressibleVoF.C
+
+LIB = $(FOAM_LIBBIN)/libincompressibleVoF

applications/solvers/modules/incompressibleVoF/Make/options

@@ -1,4 +1,6 @@
 EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
     -IincompressibleInterPhaseTransportModel/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/VoF \
@@ -6,7 +8,6 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/phaseIncompressible/lnInclude \
@@ -14,11 +15,12 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
-EXE_LIBS = \
+LIB_LIBS = \
+    -lVoFSolver \
     -lincompressibleInterPhaseTransportModel \
     -lphysicalProperties \
     -linterfaceCompression \
-    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lincompressibleTwoPhaseMixture \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lphaseIncompressibleMomentumTransportModels \

applications/solvers/modules/incompressibleVoF/alphaSuSp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,20 +23,26 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "immiscibleIncompressibleThreePhaseMixture.H"
+#include "incompressibleVoF.H"
+#include "fvcDiv.H"
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::immiscibleIncompressibleThreePhaseMixture::
-immiscibleIncompressibleThreePhaseMixture
+void Foam::solvers::incompressibleVoF::alphaSuSp
 (
-    const volVectorField& U
+    tmp<volScalarField::Internal>& Su,
+    tmp<volScalarField::Internal>& Sp
 )
-:
-    incompressibleThreePhaseMixture(U.mesh()),
-    threePhaseInterfaceProperties(*this, alpha1(), alpha2(), alpha3(), U)
-{}
+{
+    if (divergent())
+    {
+        // Phase change alpha1 source
+        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
+
+        Su = alphaSup.Su();
+        Sp = alphaSup.Sp();
+    }
+}
 
 
 // ************************************************************************* //

applications/solvers/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/Make/options

@@ -3,7 +3,6 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -11,7 +10,7 @@ EXE_INC = \
 
 LIB_LIBS = \
     -lphysicalProperties \
-    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lincompressibleTwoPhaseMixture \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lfiniteVolume \

applications/solvers/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.C

@@ -24,8 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "VoFTurbulenceDamping.H"
-#include "immiscibleIncompressibleTwoPhaseMixture.H"
-#include "interfaceProperties.H"
+#include "incompressibleTwoPhaseMixture.H"
 #include "incompressibleMomentumTransportModel.H"
 #include "fvMatrix.H"
 #include "addToRunTimeSelectionTable.H"
@@ -63,12 +62,11 @@ Foam::fv::VoFTurbulenceDamping::VoFTurbulenceDamping
     delta_("delta", dimLength, dict),
     mixture_
     (
-        mesh.lookupObject<immiscibleIncompressibleTwoPhaseMixture>
+        mesh.lookupObject<incompressibleTwoPhaseMixture>
         (
             "phaseProperties"
         )
     ),
-    interface_(mixture_),
     turbulence_
     (
         mesh.lookupType<incompressibleMomentumTransportModel>(phaseName_)
@@ -137,11 +135,13 @@ void Foam::fv::VoFTurbulenceDamping::addSup
 
     if (fieldName == "epsilon")
     {
-        eqn += interface_.fraction()*C2_*aSqrnu*turbulence_.k()()/pow4(delta_);
+        eqn += mixture_.interfaceFraction()*C2_*aSqrnu*turbulence_.k()()
+           /pow4(delta_);
     }
     else if (fieldName == "omega")
     {
-        eqn += interface_.fraction()*beta_*aSqrnu/(sqr(betaStar_)*pow4(delta_));
+        eqn += mixture_.interfaceFraction()*beta_*aSqrnu/(sqr(betaStar_)
+           *pow4(delta_));
     }
     else
     {
@@ -184,12 +184,12 @@ void Foam::fv::VoFTurbulenceDamping::addSup
 
     if (fieldName == IOobject::groupName("epsilon", phaseName_))
     {
-        eqn += interface_.fraction()*C2_*taSqrnu*turbulence_.k()()
+        eqn += mixture_.interfaceFraction()*C2_*taSqrnu*turbulence_.k()()
            /pow4(delta_);
     }
     else if (fieldName == IOobject::groupName("omega", phaseName_))
     {
-        eqn += interface_.fraction()*beta_*taSqrnu
+        eqn += mixture_.interfaceFraction()*beta_*taSqrnu
            /(sqr(betaStar_)*pow4(delta_));
     }
     else

applications/solvers/modules/incompressibleVoF/fvModels/VoFTurbulenceDamping/VoFTurbulenceDamping.H

@@ -81,8 +81,7 @@ SourceFiles
 namespace Foam
 {
 
-class immiscibleIncompressibleTwoPhaseMixture;
-class interfaceProperties;
+class incompressibleTwoPhaseMixture;
 class incompressibleMomentumTransportModel;
 
 namespace fv
@@ -109,10 +108,7 @@ class VoFTurbulenceDamping
         dimensionedScalar delta_;
 
         //- Reference to the mixture properties
-        const immiscibleIncompressibleTwoPhaseMixture& mixture_;
-
-        //- Reference to the interface properties
-        const interfaceProperties& interface_;
+        const incompressibleTwoPhaseMixture& mixture_;
 
         //- Reference to the mixture turbulence model
         const incompressibleMomentumTransportModel& turbulence_;

applications/solvers/modules/incompressibleVoF/incompressibleVoF.C

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleVoF.H"
+#include "localEulerDdtScheme.H"
+#include "CorrectPhi.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(incompressibleVoF, 0);
+    addToRunTimeSelectionTable(solver, incompressibleVoF, fvMesh);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::incompressibleVoF::incompressibleVoF(fvMesh& mesh)
+:
+    VoFSolver
+    (
+        mesh,
+        autoPtr<twoPhaseMixture>(new incompressibleTwoPhaseMixture(mesh))
+    ),
+
+    mixture(refCast<incompressibleTwoPhaseMixture>(VoFSolver::mixture)),
+
+    p
+    (
+        IOobject
+        (
+            "p",
+            runTime.name(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p_rgh + rho*buoyancy.gh
+    ),
+
+    pressureReference_
+    (
+        p,
+        p_rgh,
+        pimple.dict(),
+        false
+    ),
+
+    momentumTransport
+    (
+        U,
+        phi,
+        alphaPhi1,
+        mixture
+    )
+{
+    // Read the controls
+    read();
+
+    if (!runTime.restart() || !divergent())
+    {
+        if (correctPhi)
+        {
+            rAU = new volScalarField
+            (
+                IOobject
+                (
+                    "rAU",
+                    runTime.name(),
+                    mesh,
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh,
+                dimensionedScalar(dimTime/dimDensity, 1)
+            );
+
+            correctUphiBCs(U, phi, true);
+
+            CorrectPhi
+            (
+                phi,
+                U,
+                p_rgh,
+                surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+                geometricZeroField(),
+                pressureReference(),
+                pimple
+            );
+        }
+        else
+        {
+            correctUphiBCs(U, phi, true);
+
+            CorrectPhi
+            (
+                phi,
+                U,
+                p_rgh,
+                dimensionedScalar(dimTime/rho.dimensions(), 1),
+                geometricZeroField(),
+                pressureReference(),
+                pimple
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::incompressibleVoF::~incompressibleVoF()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleVoF::prePredictor()
+{
+    VoFSolver::prePredictor();
+
+    const dimensionedScalar& rho1 = mixture.rho1();
+    const dimensionedScalar& rho2 = mixture.rho2();
+
+    // Calculate the mass-flux from the accumulated alphaPhi1
+    rhoPhi = (alphaPhi1*(rho1 - rho2) + phi*rho2);
+}
+
+
+void Foam::solvers::incompressibleVoF::thermophysicalPredictor()
+{}
+
+
+// ************************************************************************* //

applications/solvers/modules/incompressibleVoF/incompressibleVoF.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+Class
+    Foam::solvers::incompressibleVoF
+
+Description
+    Solver module for 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.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Either mixture or two-phase transport modelling may be selected.  In the
+    mixture approach a single laminar, RAS or LES model is selected to model the
+    momentum stress.  In the Euler-Euler two-phase approach separate laminar,
+    RAS or LES selected models are selected for each of the phases.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient and steady simulations.
+
+    Optional fvModels and fvConstraints are provided to enhance the simulation
+    in many ways including adding various sources, Lagrangian
+    particles, surface film etc. and constraining or limiting the solution.
+
+SourceFiles
+    incompressibleVoF.C
+
+See also
+    Foam::solvers::VoFSolver
+    Foam::solvers::compressibleVoF
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef incompressibleVoF_H
+#define incompressibleVoF_H
+
+#include "VoFSolver.H"
+#include "incompressibleTwoPhaseMixture.H"
+#include "incompressibleInterPhaseTransportModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class incompressibleVoF Declaration
+\*---------------------------------------------------------------------------*/
+
+class incompressibleVoF
+:
+    public VoFSolver
+{
+
+protected:
+
+    // Phase properties
+
+        //- The compressible two-phase mixture
+        incompressibleTwoPhaseMixture& mixture;
+
+
+    // Thermophysical properties
+
+        //- Static pressure field
+        volScalarField p;
+
+
+    // Pressure reference
+
+        //- Pressure reference
+        Foam::pressureReference pressureReference_;
+
+
+    // Momentum transport
+
+        //- Momentum transport model
+        incompressibleInterPhaseTransportModel momentumTransport;
+
+
+    // Protected Member Functions
+
+        //- Return the pressure reference
+        virtual const Foam::pressureReference& pressureReference() const
+        {
+            return pressureReference_;
+        }
+
+        //- Is the flow divergent?
+        //  i.e. includes phase-fraction sources
+        virtual bool divergent()
+        {
+            return fvModels().addsSupToField(alpha1.name());
+        }
+
+        //- Calculate the alpha equation sources
+        virtual void alphaSuSp
+        (
+            tmp<volScalarField::Internal>& Su,
+            tmp<volScalarField::Internal>& Sp
+        );
+
+        //- Return the momentum equation stress term
+        virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U)
+        {
+            return momentumTransport.divDevTau(rho, U);
+        }
+
+
+public:
+
+    //- Runtime type information
+    TypeName("incompressibleVoF");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        incompressibleVoF(fvMesh& mesh);
+
+        //- Disallow default bitwise copy construction
+        incompressibleVoF(const incompressibleVoF&) = delete;
+
+
+    //- Destructor
+    virtual ~incompressibleVoF();
+
+
+    // Member Functions
+
+        //- Called at the start of the PIMPLE loop
+        virtual void prePredictor();
+
+        //- Construct and solve the energy equation,
+        //  convert to temperature
+        //  and update thermophysical and transport properties
+        virtual void thermophysicalPredictor();
+
+        //- Construct and solve the pressure equation in the PISO loop
+        virtual void pressureCorrector();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const incompressibleVoF&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/incompressibleVoF/pressureCorrector.C

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleVoF.H"
+#include "constrainHbyA.H"
+#include "constrainPressure.H"
+#include "adjustPhi.H"
+#include "findRefCell.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleVoF::pressureCorrector()
+{
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    if (rAU.valid())
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+
+    while (pimple.correct())
+    {
+        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
+        (
+            (
+                interface.surfaceTensionForce()
+              - buoyancy.ghf*fvc::snGrad(rho)
+            )*rAUf*mesh.magSf()
+        );
+
+        phiHbyA += phig;
+
+        // Update the pressure BCs to ensure flux consistency
+        constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
+
+        // Cache the phase change pressure source
+        fvScalarMatrix Sp_rgh
+        (
+            fvModels().source
+            (
+                volScalarField::New
+                (
+                    "1",
+                    mesh,
+                    dimensionedScalar(dimless/dimPressure, 1)
+                ),
+                p_rgh
+            )
+        );
+
+        while (pimple.correctNonOrthogonal())
+        {
+            fvScalarMatrix p_rghEqn
+            (
+                fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
+             == Sp_rgh
+            );
+
+            p_rghEqn.setReference
+            (
+                pressureReference().refCell(),
+                getRefCellValue(p_rgh, pressureReference().refCell())
+            );
+
+            p_rghEqn.solve();
+
+            if (pimple.finalNonOrthogonalIter())
+            {
+                phi = phiHbyA + p_rghEqn.flux();
+
+                p_rgh.relax();
+
+                U = HbyA
+                  + rAU()*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+                U.correctBoundaryConditions();
+                fvConstraints().constrain(U);
+            }
+        }
+
+        continuityErrors();
+
+        // Correct Uf if the mesh is moving
+        fvc::correctUf(Uf, U, phi, MRF);
+
+        // Make the fluxes relative to the mesh motion
+        fvc::makeRelative(phi, U);
+
+        p == p_rgh + rho*buoyancy.gh;
+
+        if (p_rgh.needReference())
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pressureReference().refValue()
+              - getRefCellValue(p, pressureReference().refCell())
+            );
+            p_rgh = p - rho*buoyancy.gh;
+        }
+    }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
+
+    tUEqn.clear();
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/cavitatingFoam/CourantNo.H

@@ -44,7 +44,7 @@ scalar acousticCoNum = 0;
     (
         fvc::surfaceSum
         (
-            fvc::interpolate(scalar(1)/sqrt(psi))*mesh.magSf()
+            fvc::interpolate(scalar(1)/sqrt(mixture.psi()))*mesh.magSf()
         )().primitiveField()/mesh.V().field()
     )*runTime.deltaTValue();
 

applications/solvers/multiphase/cavitatingFoam/Make/files

@@ -1,3 +1,4 @@
+cavitatingTwoPhaseMixture/cavitatingTwoPhaseMixture.C
 cavitatingFoam.C
 
 EXE = $(FOAM_APPBIN)/cavitatingFoam

applications/solvers/multiphase/cavitatingFoam/Make/options

@@ -1,23 +1,20 @@
 EXE_INC = \
+    -IcavitatingTwoPhaseMixture \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/VoF \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceCompression/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lphysicalProperties \
-    -linterfaceCompression \
-    -lincompressibleTwoPhaseMixture \
+    -ltwoPhaseMixture \
+    -lbarotropicCompressibilityModel \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
-    -lbarotropicCompressibilityModel \
     -lfiniteVolume \
     -lmeshTools \
     -lfvModels \

applications/solvers/multiphase/cavitatingFoam/alphavPsi.H

@@ -1,20 +0,0 @@
-{
-    alphav =
-        max
-        (
-            min
-            (
-                (rho - rholSat)/(rhovSat - rholSat),
-                scalar(1)
-            ),
-            scalar(0)
-        );
-    alphal = 1.0 - alphav;
-
-    Info<< "max-min alphav: " << max(alphav).value()
-        << " " << min(alphav).value() << endl;
-
-    psiModel->correct();
-
-    // Info<< "min a: " << 1.0/sqrt(max(psi)).value() << endl;
-}

applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C

@@ -34,8 +34,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "barotropicCompressibilityModel.H"
-#include "incompressibleTwoPhaseMixture.H"
+#include "cavitatingTwoPhaseMixture.H"
 #include "incompressibleMomentumTransportModels.H"
 #include "CorrectPhi.H"
 #include "pimpleControl.H"
@@ -90,7 +89,8 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "rhoEqn.H"
-            #include "alphavPsi.H"
+
+            mixture.correct();
 
             if (pimple.predictTransport())
             {

applications/solvers/multiphase/cavitatingFoam/cavitatingTwoPhaseMixture/cavitatingTwoPhaseMixture.C

@@ -0,0 +1,247 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cavitatingTwoPhaseMixture.H"
+#include "barotropicCompressibilityModel.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cavitatingTwoPhaseMixture, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cavitatingTwoPhaseMixture::cavitatingTwoPhaseMixture
+(
+    const fvMesh& mesh
+)
+:
+    twoPhaseMixture(mesh),
+
+    alphav_(alpha1()),
+    alphal_(alpha2()),
+
+    nuModelv_(viscosityModel::New(mesh, phase1Name())),
+    nuModell_(viscosityModel::New(mesh, phase2Name())),
+
+    rhov_("rho", dimDensity, nuModelv_()),
+    rhol_("rho", dimDensity, nuModell_()),
+
+    nu_
+    (
+        IOobject
+        (
+            "nu",
+            mesh.time().name(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar(dimViscosity, 0),
+        calculatedFvPatchScalarField::typeName
+    ),
+
+    thermodynamicProperties_
+    (
+        IOobject
+        (
+            "thermodynamicProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+
+    psil_
+    (
+        "psil",
+        dimCompressibility,
+        thermodynamicProperties_
+    ),
+
+    rholSat_
+    (
+        "rholSat",
+        dimDensity,
+        thermodynamicProperties_
+    ),
+
+    psiv_
+    (
+        "psiv",
+        dimCompressibility,
+        thermodynamicProperties_
+    ),
+
+    pSat_
+    (
+        "pSat",
+        dimPressure,
+        thermodynamicProperties_
+    ),
+
+    rhovSat_("rhovSat", psiv_*pSat_),
+
+    rhol0_("rhol0", rholSat_ - pSat_*psil_),
+
+    rhoMin_
+    (
+        "rhoMin",
+        dimDensity,
+        thermodynamicProperties_
+    ),
+
+    p_
+    (
+        IOobject
+        (
+            "p",
+            mesh.time().name(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    ),
+
+    psiModel_
+    (
+        barotropicCompressibilityModel::New
+        (
+            thermodynamicProperties_,
+            alphav_
+        )
+    ),
+
+    psi_(psiModel_->psi()),
+
+    rho_
+    (
+        IOobject
+        (
+            "rho",
+            mesh.time().name(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    )
+{
+    alphav_.oldTime();
+
+    rho_ = max
+    (
+        psi_*p_
+      + alphal_*rhol0_
+      + ((alphav_*psiv_ + alphal_*psil_) - psi_)*pSat_,
+        rhoMin_
+    );
+
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cavitatingTwoPhaseMixture::~cavitatingTwoPhaseMixture()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::cavitatingTwoPhaseMixture::read()
+{
+    if (twoPhaseMixture::read())
+    {
+        nuModelv_->lookup("rho") >> rhov_;
+        nuModell_->lookup("rho") >> rhol_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::cavitatingTwoPhaseMixture::correctPressure()
+{
+    p_ =
+    (
+        rho_
+      - alphal_*rhol0_
+      - ((alphav_*psiv_ + alphal_*psil_) - psi_)*pSat_
+    )/psi_;
+
+    p_.correctBoundaryConditions();
+}
+
+
+void Foam::cavitatingTwoPhaseMixture::correct()
+{
+    rho_ == max(rho_, rhoMin_);
+
+    alphav_ =
+        max
+        (
+            min
+            (
+                (rho_ - rholSat_)/(rhovSat_ - rholSat_),
+                scalar(1)
+            ),
+            scalar(0)
+        );
+    alphal_ = 1.0 - alphav_;
+
+    Info<< "max-min alphav: " << max(alphav_).value()
+        << " " << min(alphav_).value() << endl;
+
+    psiModel_->correct();
+
+    nuModelv_->correct();
+    nuModell_->correct();
+
+    const volScalarField limitedAlphav
+    (
+        "limitedAlphav",
+        min(max(alphav_, scalar(0)), scalar(1))
+    );
+
+    // Mixture kinematic viscosity calculated from mixture dynamic viscosity
+    nu_ =
+    (
+        limitedAlphav*rhov_*nuModelv_->nu()
+      + (scalar(1) - limitedAlphav)*rhol_*nuModell_->nu()
+    )/(limitedAlphav*rhov_ + (scalar(1) - limitedAlphav)*rhol_);
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/cavitatingFoam/cavitatingTwoPhaseMixture/cavitatingTwoPhaseMixture.H

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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/>.
+
+Class
+    Foam::cavitatingTwoPhaseMixture
+
+Description
+    Class to represent a mixture of two constant density phases
+
+SourceFiles
+    cavitatingTwoPhaseMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cavitatingTwoPhaseMixture_H
+#define cavitatingTwoPhaseMixture_H
+
+#include "twoPhaseMixture.H"
+#include "viscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class barotropicCompressibilityModel;
+
+/*---------------------------------------------------------------------------*\
+                Class cavitatingTwoPhaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class cavitatingTwoPhaseMixture
+:
+    public twoPhaseMixture,
+    public viscosity
+{
+    // Private Data
+
+        //- Vapour phase fraction reference (alpha1 from twoPhaseMixture)
+        volScalarField& alphav_;
+
+        //- liquid phase fraction reference (alpha2 from twoPhaseMixture)
+        volScalarField& alphal_;
+
+        //- Vapour viscosity model
+        autoPtr<viscosityModel> nuModelv_;
+
+        //- Liquid viscosity model
+        autoPtr<viscosityModel> nuModell_;
+
+        //- Constant vapour density
+        dimensionedScalar rhov_;
+
+        //- Constant liquid density
+        dimensionedScalar rhol_;
+
+        //- Mixture viscosity
+        volScalarField nu_;
+
+        //- Thermophysical properties dictionary
+        IOdictionary thermodynamicProperties_;
+
+        //- Liquid compressibility
+        dimensionedScalar psil_;
+
+        //- Saturation liquid density
+        dimensionedScalar rholSat_;
+
+        //- Vapour compressibility
+        dimensionedScalar psiv_;
+
+        //- Saturation pressure
+        dimensionedScalar pSat_;
+
+        //- Saturation vapour density
+        dimensionedScalar rhovSat_;
+
+        //- rholSat_ - pSat_*psil_
+        dimensionedScalar rhol0_;
+
+        //- Minimum density for numerical stability
+        dimensionedScalar rhoMin_;
+
+        //- Static pressure field
+        volScalarField p_;
+
+        autoPtr<barotropicCompressibilityModel> psiModel_;
+
+        //- Reference to the mixture compressibility field
+        const volScalarField& psi_;
+
+        //- Mixture density field
+        volScalarField rho_;
+
+
+public:
+
+    TypeName("cavitatingTwoPhaseMixture");
+
+
+    // Constructors
+
+        //- Construct from a mesh
+        cavitatingTwoPhaseMixture(const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~cavitatingTwoPhaseMixture();
+
+
+    // Member Functions
+
+        //- Return the static pressure field
+        volScalarField& p()
+        {
+            return p_;
+        }
+
+        //- Return the mixture density field
+        virtual const volScalarField& rho() const
+        {
+            return rho_;
+        }
+
+        //- Return the mixture density field
+        volScalarField& rho()
+        {
+            return rho_;
+        }
+
+        //- Return the mixture compressibility
+        const volScalarField& psi() const
+        {
+            return psi_;
+        }
+
+        //- Return the kinematic laminar viscosity
+        virtual tmp<volScalarField> nu() const
+        {
+            return nu_;
+        }
+
+        //- Return the laminar viscosity for patch
+        virtual tmp<scalarField> nu(const label patchi) const
+        {
+            return nu_.boundaryField()[patchi];
+        }
+
+        //- Correct the pressure from mixture density
+        void correctPressure();
+
+        //- Correct the mixture density and laminar viscosity
+        virtual void correct();
+
+        //- Read base phaseProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/cavitatingFoam/continuityErrs.H

@@ -1,22 +0,0 @@
-{
-    volScalarField thermoRho = psi*p + alphal*rhol0;
-
-    dimensionedScalar totalMass = fvc::domainIntegrate(rho);
-
-    scalar sumLocalContErr =
-    (
-        fvc::domainIntegrate(mag(rho - thermoRho))/totalMass
-    ).value();
-
-    scalar globalContErr =
-    (
-        fvc::domainIntegrate(rho - thermoRho)/totalMass
-    ).value();
-
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr
-        << endl;
-}

applications/solvers/multiphase/cavitatingFoam/createFields.H

@@ -1,31 +1,9 @@
-#include "readThermodynamicProperties.H"
+Info<< "Reading phaseProperties\n" << endl;
 
-Info<< "Reading field p\n" << endl;
-volScalarField p
-(
-    IOobject
-    (
-        "p",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
+cavitatingTwoPhaseMixture mixture(mesh);
 
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
+volScalarField& p(mixture.p());
+volScalarField& rho(mixture.rho());
 
 Info<< "Reading field U\n" << endl;
 volVectorField U
@@ -59,33 +37,6 @@ surfaceScalarField rhoPhi
     fvc::interpolate(rho)*phi
 );
 
-Info<< "Reading phaseProperties\n" << endl;
-
-incompressibleTwoPhaseMixture mixture(mesh);
-
-volScalarField& alphav(mixture.alpha1());
-alphav.oldTime();
-
-volScalarField& alphal(mixture.alpha2());
-
-Info<< "Creating compressibilityModel\n" << endl;
-autoPtr<barotropicCompressibilityModel> psiModel =
-    barotropicCompressibilityModel::New
-    (
-        thermodynamicProperties,
-        alphav
-    );
-
-const volScalarField& psi = psiModel->psi();
-
-rho == max
-(
-    psi*p
-  + alphal*rhol0
-  + ((alphav*psiv + alphal*psil) - psi)*pSat,
-    rhoMin
-);
-
 mesh.schemes().setFluxRequired(p.name());
 mesh.schemes().setFluxRequired(rho.name());
 

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -1,12 +1,7 @@
 {
     if (pimple.nCorr() == 1)
     {
-        p =
-        (
-            rho
-          - alphal*rhol0
-          - ((alphav*psiv + alphal*psil) - psi)*pSat
-        )/psi;
+        mixture.correctPressure();
     }
 
     surfaceScalarField rhof("rhof", fvc::interpolate(rho));
@@ -22,14 +17,14 @@
 
     phi -= phiGradp/rhof;
 
-    volScalarField rho0(rho - psi*p);
+    const volScalarField rho0(rho - mixture.psi()*p);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvc::ddt(rho)
-          + psi*correction(fvm::ddt(p))
+          + mixture.psi()*correction(fvm::ddt(p))
           + fvc::div(phi, rho)
           + fvc::div(phiGradp)
           - fvm::laplacian(rhorAUf, p)
@@ -46,18 +41,11 @@
     Info<< "Predicted p max-min : " << max(p).value()
         << " " << min(p).value() << endl;
 
-    rho == max(rho0 + psi*p, rhoMin);
+    rho == rho0 + mixture.psi()*p;
 
-    #include "alphavPsi.H"
+    mixture.correct();
 
-    p =
-    (
-        rho
-      - alphal*rhol0
-      - ((alphav*psiv + alphal*psil) - psi)*pSat
-    )/psi;
-
-    p.correctBoundaryConditions();
+    mixture.correctPressure();
 
     Info<< "Phase-change corrected p max-min : " << max(p).value()
         << " " << min(p).value() << endl;

applications/solvers/multiphase/cavitatingFoam/readThermodynamicProperties.H

@@ -1,52 +0,0 @@
-Info<< "Reading thermodynamicProperties\n" << endl;
-
-IOdictionary thermodynamicProperties
-(
-    IOobject
-    (
-        "thermodynamicProperties",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ_IF_MODIFIED,
-        IOobject::NO_WRITE
-    )
-);
-
-dimensionedScalar psil
-(
-    "psil",
-    dimCompressibility,
-    thermodynamicProperties
-);
-
-dimensionedScalar rholSat
-(
-    "rholSat",
-    dimDensity,
-    thermodynamicProperties
-);
-
-dimensionedScalar psiv
-(
-    "psiv",
-    dimCompressibility,
-    thermodynamicProperties
-);
-
-dimensionedScalar pSat
-(
-    "pSat",
-    dimPressure,
-    thermodynamicProperties
-);
-
-dimensionedScalar rhovSat("rhovSat", psiv*pSat);
-
-dimensionedScalar rhol0("rhol0", rholSat - pSat*psil);
-
-dimensionedScalar rhoMin
-(
-    "rhoMin",
-    dimDensity,
-    thermodynamicProperties
-);

applications/solvers/multiphase/cavitatingFoam/rhoEqn.H

@@ -11,6 +11,4 @@
 
     Info<< "max-min rho: " << max(rho).value()
         << " " << min(rho).value() << endl;
-
-    rho == max(rho, rhoMin);
 }

applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H

@@ -68,6 +68,7 @@
             << endl;
     }
 
+    mixture.correct();
+
     rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
-    rho = mixture.rho();
 }

applications/solvers/multiphase/driftFluxFoam/createFields.H

@@ -38,20 +38,7 @@ Info<< "Reading incompressibleTwoPhaseInteractingMixture\n" << endl;
 incompressibleTwoPhaseInteractingMixture mixture(U, phi, g);
 
 volScalarField& alpha1(mixture.alpha1());
-
-// Mixture density
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.name(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    mixture.rho()
-);
+const volScalarField& rho(mixture.rho());
 
 // Mass flux
 surfaceScalarField rhoPhi

applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.C

@@ -64,7 +64,19 @@ incompressibleTwoPhaseInteractingMixture
     ),
     alphaMax_(lookupOrDefault("alphaMax", 1.0)),
 
-    MRF_(U.mesh()),
+    rho_
+    (
+        IOobject
+        (
+            "rho",
+            U_.time().name(),
+            U_.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        U_.mesh(),
+        dimensionedScalar("rho", dimDensity, 0)
+    ),
 
     mu_
     (
@@ -79,6 +91,8 @@ incompressibleTwoPhaseInteractingMixture
         calculatedFvPatchScalarField::typeName
     ),
 
+    MRF_(U.mesh()),
+
     UdmModel_(relativeVelocityModel::New(*this, *this, g))
 {
     correct();
@@ -156,11 +170,10 @@ Foam::incompressibleTwoPhaseInteractingMixture::U() const
     return U_;
 }
 
-
-const Foam::IOMRFZoneList&
-Foam::incompressibleTwoPhaseInteractingMixture::MRF() const
+const Foam::volScalarField&
+Foam::incompressibleTwoPhaseInteractingMixture::rho() const
 {
-    return MRF_;
+    return rho_;
 }
 
 
@@ -178,33 +191,24 @@ Foam::incompressibleTwoPhaseInteractingMixture::mu(const label patchi) const
 }
 
 
-Foam::tmp<Foam::volScalarField>
-Foam::incompressibleTwoPhaseInteractingMixture::rho() const
-{
-    return alpha1()*rhod_ + alpha2()*rhoc_;
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::incompressibleTwoPhaseInteractingMixture::rho(const label patchi) const
-{
-    return
-        alpha1().boundaryField()[patchi]*rhod_.value()
-      + alpha2().boundaryField()[patchi]*rhoc_.value();
-}
-
-
 Foam::tmp<Foam::volScalarField>
 Foam::incompressibleTwoPhaseInteractingMixture::nu() const
 {
-    return mu_/rho();
+    return mu_/rho_;
 }
 
 
 Foam::tmp<Foam::scalarField>
 Foam::incompressibleTwoPhaseInteractingMixture::nu(const label patchi) const
 {
-    return mu_.boundaryField()[patchi]/rho(patchi);
+    return mu_.boundaryField()[patchi]/rho_.boundaryField()[patchi];
+}
+
+
+const Foam::IOMRFZoneList&
+Foam::incompressibleTwoPhaseInteractingMixture::MRF() const
+{
+    return MRF_;
 }
 
 
@@ -224,6 +228,7 @@ Foam::incompressibleTwoPhaseInteractingMixture::divTauDm() const
 
 void Foam::incompressibleTwoPhaseInteractingMixture::correct()
 {
+    rho_ = alpha1()*rhod_ + alpha2()*rhoc_;
     mu_ = muModel_->mu(rhoc_*nucModel_->nu(), U_);
     UdmModel_->correct();
 }

applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.H

@@ -75,11 +75,13 @@ class incompressibleTwoPhaseInteractingMixture
         //- Optional maximum dispersed phase-fraction (e.g. packing limit)
         scalar alphaMax_;
 
-        //- Optional MRF zones
-        IOMRFZoneList MRF_;
+        volScalarField rho_;
 
         volScalarField mu_;
 
+        //- Optional MRF zones
+        IOMRFZoneList MRF_;
+
         //- Dispersed phase relative velocity model
         autoPtr<relativeVelocityModel> UdmModel_;
 
@@ -134,8 +136,8 @@ public:
         //- Return const-access to the mixture velocity
         const volVectorField& U() const;
 
-        //- Return MRF zones
-        const IOMRFZoneList& MRF() const;
+        //- Return the mixture density
+        virtual const volScalarField& rho() const;
 
         //- Return the dynamic mixture viscosity
         tmp<volScalarField> mu() const;
@@ -143,18 +145,15 @@ public:
         //- Return the dynamic mixture viscosity for patch
         virtual tmp<scalarField> mu(const label patchi) const;
 
-        //- Return the mixture density
-        virtual tmp<volScalarField> rho() const;
-
-        //- Return the mixture density for patch
-        virtual tmp<scalarField> rho(const label patchi) const;
-
         //- Return the mixture viscosity
         virtual tmp<volScalarField> nu() const;
 
         //- Return the mixture viscosity for patch
         virtual tmp<scalarField> nu(const label patchi) const;
 
+        //- Return MRF zones
+        const IOMRFZoneList& MRF() const;
+
         //- Return the diffusion velocity of the dispersed phase
         const volVectorField& Udm() const;
 

applications/solvers/multiphase/interFoam/Make/files

@@ -1,3 +0,0 @@
-interFoam.C
-
-EXE = $(FOAM_APPBIN)/interFoam

applications/solvers/multiphase/interFoam/UEqn.H

@@ -1,32 +0,0 @@
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      + MRF.DDt(rho, U)
-      + turbulence.divDevTau(rho, U)
-     ==
-        // phaseChange.SU(rho, rhoPhi, U)
-       fvModels.source(rho, U)
-    );
-
-    UEqn.relax();
-
-    fvConstraints.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                    mixture.surfaceTensionForce()
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                ) * mesh.magSf()
-            )
-        );
-
-        fvConstraints.constrain(U);
-    }

applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H

@@ -1,50 +0,0 @@
-if (nAlphaSubCycles > 1)
-{
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-    tmp<volScalarField> trSubDeltaT;
-
-    if (LTS)
-    {
-        trSubDeltaT =
-            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
-    }
-
-    // Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1
-    tmp<surfaceScalarField> talphaPhi1
-    (
-        surfaceScalarField::New
-        (
-            "alphaPhi1",
-            mesh,
-            dimensionedScalar(alphaPhi1.dimensions(), 0)
-        )
-    );
-
-    List<volScalarField*> alphaPtrs({&alpha1, &alpha2});
-
-    for
-    (
-        subCycle<volScalarField, subCycleFields> alphaSubCycle
-        (
-            alphaPtrs,
-            nAlphaSubCycles
-        );
-        !(++alphaSubCycle).end();
-    )
-    {
-        #include "alphaEqn.H"
-        talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
-    }
-
-    alphaPhi1 = talphaPhi1();
-}
-else
-{
-    #include "alphaEqn.H"
-}
-
-mixture.correct();
-rho == mixture.rho();
-
-// Calculate the mass-flux from the accumulated alphaPhi1
-rhoPhi = (alphaPhi1*(rho1 - rho2) + phi*rho2);

applications/solvers/multiphase/interFoam/correctPhi.H

@@ -1,37 +0,0 @@
-// Calculate absolute flux
-// from the mapped surface velocity
-phi = mesh.Sf() & Uf();
-
-correctUphiBCs(U, phi, true);
-
-if (divU.valid())
-{
-    CorrectPhi
-    (
-        phi,
-        U,
-        p_rgh,
-        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-        divU(),
-        pressureReference,
-        pimple
-    );
-}
-else
-{
-    CorrectPhi
-    (
-        phi,
-        U,
-        p_rgh,
-        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-        geometricZeroField(),
-        pressureReference,
-        pimple
-    );
-}
-
-#include "continuityErrs.H"
-
-// Make the flux relative to the mesh motion
-fvc::makeRelative(phi, U);

applications/solvers/multiphase/interFoam/createFieldRefs.H

@@ -1,7 +0,0 @@
-volScalarField& alpha2(mixture.alpha2());
-
-const dimensionedScalar& rho1 = mixture.rho1();
-const dimensionedScalar& rho2 = mixture.rho2();
-
-//twoPhaseChangeModel& phaseChange = phaseChangePtr();
-tmp<volScalarField> rAU;

applications/solvers/multiphase/interFoam/createFields.H

@@ -1,105 +0,0 @@
-#include "createRDeltaT.H"
-
-Info<< "Reading field p_rgh\n" << endl;
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "createPhi.H"
-
-
-Info<< "Reading phaseProperties\n" << endl;
-immiscibleIncompressibleTwoPhaseMixture mixture(U);
-
-
-// Need to store rho for ddt(rho, U)
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.name(),
-        mesh,
-        IOobject::READ_IF_PRESENT
-    ),
-    mixture.rho()
-);
-rho.oldTime();
-
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-volScalarField p
-(
-    IOobject
-    (
-        "p",
-        runTime.name(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    p_rgh + rho*gh
-);
-
-pressureReference pressureReference(p, p_rgh, pimple.dict());
-
-if (p_rgh.needReference())
-{
-    p += dimensionedScalar
-    (
-        "p",
-        p.dimensions(),
-        pressureReference.refValue()
-      - getRefCellValue(p, pressureReference.refCell())
-    );
-    p_rgh = p - rho*gh;
-}
-
-mesh.schemes().setFluxRequired(p_rgh.name());
-
-volScalarField& alpha1(mixture.alpha1());
-mesh.schemes().setFluxRequired(alpha1.name());
-
-
-#include "createAlphaFluxes.H"
-
-// Construct incompressible turbulence model
-incompressibleInterPhaseTransportModel turbulence
-(
-    U,
-    phi,
-    alphaPhi1,
-    mixture
-);
-
-
-#include "createMRF.H"
-#include "createFvModels.H"
-#include "createFvConstraints.H"

applications/solvers/multiphase/interFoam/initCorrectPhi.H

@@ -1,53 +0,0 @@
-if
-(
-    !runTime.restart()
- || !fvModels.addsSupToField(alpha1.name())
-)
-{
-    if (correctPhi)
-    {
-        rAU = new volScalarField
-        (
-            IOobject
-            (
-                "rAU",
-                runTime.name(),
-                mesh,
-                IOobject::READ_IF_PRESENT,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            dimensionedScalar(dimTime/dimDensity, 1)
-        );
-
-        correctUphiBCs(U, phi, true);
-
-        CorrectPhi
-        (
-            phi,
-            U,
-            p_rgh,
-            surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-            geometricZeroField(),
-            pressureReference,
-            pimple
-        );
-    }
-    else
-    {
-        correctUphiBCs(U, phi, true);
-
-        CorrectPhi
-        (
-            phi,
-            U,
-            p_rgh,
-            dimensionedScalar(dimTime/rho.dimensions(), 1),
-            geometricZeroField(),
-            pressureReference,
-            pimple
-        );
-    }
-}
-
-#include "continuityErrs.H"

applications/solvers/multiphase/interFoam/interFoam.C

@@ -1,226 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 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,
-    with optional mesh motion and mesh topology changes including adaptive
-    re-meshing.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "interfaceCompression.H"
-#include "CMULES.H"
-#include "EulerDdtScheme.H"
-#include "localEulerDdtScheme.H"
-#include "CrankNicolsonDdtScheme.H"
-#include "subCycle.H"
-#include "immiscibleIncompressibleTwoPhaseMixture.H"
-#include "incompressibleInterPhaseTransportModel.H"
-#include "pimpleControl.H"
-#include "pressureReference.H"
-#include "fvModels.H"
-#include "fvConstraints.H"
-#include "CorrectPhi.H"
-#include "fvcSmooth.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "initContinuityErrs.H"
-    #include "createDyMControls.H"
-    #include "createFields.H"
-    #include "createFieldRefs.H"
-    #include "initCorrectPhi.H"
-    #include "createUfIfPresent.H"
-
-    if (!LTS)
-    {
-        #include "CourantNo.H"
-        #include "setInitialDeltaT.H"
-    }
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (pimple.run(runTime))
-    {
-        #include "readDyMControls.H"
-
-        if (LTS)
-        {
-            #include "setRDeltaT.H"
-        }
-        else
-        {
-            #include "CourantNo.H"
-            #include "alphaCourantNo.H"
-            #include "setDeltaT.H"
-        }
-
-        fvModels.preUpdateMesh();
-
-        // Store divU from the previous mesh so that it can be mapped
-        // and used in correctPhi to ensure the corrected phi has the
-        // same divergence
-        tmp<volScalarField> divU;
-
-        if
-        (
-            correctPhi
-         && fvModels.addsSupToField(alpha1.name())
-         && mesh.topoChanged()
-        )
-        {
-            // Construct and register divU for correctPhi
-            divU = new volScalarField
-            (
-                "divU0",
-                fvc::div(fvc::absolute(phi, U))
-            );
-        }
-
-        // Update the mesh for topology change, mesh to mesh mapping
-        bool topoChanged = mesh.update();
-
-        // Do not apply previous time-step mesh compression flux
-        // if the mesh topology changed
-        if (topoChanged)
-        {
-            talphaPhi1Corr0.clear();
-        }
-
-        runTime++;
-
-        Info<< "Time = " << runTime.userTimeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            if (pimple.firstPimpleIter() || pimple.moveMeshOuterCorrectors())
-            {
-                if
-                (
-                    correctPhi
-                 && fvModels.addsSupToField(alpha1.name())
-                 && !divU.valid()
-                )
-                {
-                    // Construct and register divU for correctPhi
-                    divU = new volScalarField
-                    (
-                        "divU0",
-                        fvc::div(fvc::absolute(phi, U))
-                    );
-                }
-
-                // Move the mesh
-                mesh.move();
-
-                if (mesh.changing())
-                {
-                    gh = (g & mesh.C()) - ghRef;
-                    ghf = (g & mesh.Cf()) - ghRef;
-
-                    MRF.update();
-
-                    if (correctPhi)
-                    {
-                        #include "correctPhi.H"
-                    }
-
-                    mixture.correct();
-
-                    if (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
-                }
-
-                divU.clear();
-            }
-
-            fvModels.correct();
-
-            surfaceScalarField rhoPhi
-            (
-                IOobject
-                (
-                    "rhoPhi",
-                    runTime.name(),
-                    mesh
-                ),
-                mesh,
-                dimensionedScalar(dimMass/dimTime, 0)
-            );
-
-            #include "alphaControls.H"
-            #include "alphaEqnSubCycle.H"
-
-            turbulence.correctPhasePhi();
-
-            mixture.correct();
-
-            if (pimple.predictTransport())
-            {
-                turbulence.predict();
-            }
-
-            #include "UEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.correctTransport())
-            {
-                turbulence.correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/Make/files

@@ -1,6 +0,0 @@
-incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.C
-threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
-immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.C
-interMixingFoam.C
-
-EXE = $(FOAM_APPBIN)/interMixingFoam

applications/solvers/multiphase/interFoam/interMixingFoam/Make/options

@@ -1,28 +0,0 @@
-EXE_INC = \
-    -I. \
-    -I.. \
-    -I$(LIB_SRC)/physicalProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/VoF \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceCompression/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/twoPhaseProperties/alphaContactAngle/alphaContactAngle \
-    -IimmiscibleIncompressibleThreePhaseMixture \
-    -IincompressibleThreePhaseMixture \
-    -IthreePhaseInterfaceProperties \
-    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-    -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    -lphysicalProperties \
-    -linterfaceCompression \
-    -ltwoPhaseProperties \
-    -lmomentumTransportModels \
-    -lincompressibleMomentumTransportModels \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lfvModels \
-    -lfvConstraints \
-    -lsampling

applications/solvers/multiphase/interFoam/interMixingFoam/UEqn.H

@@ -1,31 +0,0 @@
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      + MRF.DDt(rho, U)
-      + turbulence->divDevTau(rho, U)
-     ==
-        fvModels.source(rho, U)
-    );
-
-    UEqn.relax();
-
-    fvConstraints.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                    mixture.surfaceTensionForce()
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                ) * mesh.magSf()
-            )
-        );
-
-        fvConstraints.constrain(U);
-    }

applications/solvers/multiphase/interFoam/interMixingFoam/alphaControls.H

@@ -1,5 +0,0 @@
-const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
-
-label nAlphaCorr(alphaControls.lookup<label>("nAlphaCorr"));
-
-label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H

@@ -1,211 +0,0 @@
-{
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
-
-    surfaceScalarField phir
-    (
-        IOobject
-        (
-            "phir",
-            runTime.name(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mixture.cAlpha()*mag(phi/mesh.magSf())*mixture.nHatf()
-    );
-
-    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
-    {
-        // Split the limiter into a surfaceScalarField
-        surfaceScalarField lambda
-        (
-            IOobject
-            (
-                "lambda",
-                mesh.time().name(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar(dimless, 1)
-        );
-
-        // Create the complete convection flux for alpha1
-        surfaceScalarField alphaPhi1
-        (
-            fvc::flux
-            (
-                phi,
-                alpha1,
-                alphaScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(-phir, alpha2, alpharScheme),
-                alpha1,
-                alpharScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(-phir, alpha3, alpharScheme),
-                alpha1,
-                alpharScheme
-            )
-        );
-
-        // Create the bounded (upwind) flux for alpha1
-        surfaceScalarField alphaPhi1BD
-        (
-            upwind<scalar>(mesh, phi).flux(alpha1)
-        );
-
-        // Calculate the flux correction for alpha1
-        alphaPhi1 -= alphaPhi1BD;
-
-        // Calculate the limiter for alpha1
-        if (LTS)
-        {
-            const volScalarField& rDeltaT =
-                fv::localEulerDdt::localRDeltaT(mesh);
-
-            MULES::limiter
-            (
-                lambda,
-                rDeltaT,
-                geometricOneField(),
-                alpha1,
-                alphaPhi1BD,
-                alphaPhi1,
-                zeroField(),
-                zeroField(),
-                oneField(),
-                zeroField()
-            );
-        }
-        else
-        {
-            MULES::limiter
-            (
-                lambda,
-                1.0/runTime.deltaT().value(),
-                geometricOneField(),
-                alpha1,
-                alphaPhi1BD,
-                alphaPhi1,
-                zeroField(),
-                zeroField(),
-                oneField(),
-                zeroField()
-            );
-        }
-
-        // Create the complete flux for alpha2
-        surfaceScalarField alphaPhi2
-        (
-            fvc::flux
-            (
-                phi,
-                alpha2,
-                alphaScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(phir, alpha1, alpharScheme),
-                alpha2,
-                alpharScheme
-            )
-        );
-
-        // Create the bounded (upwind) flux for alpha2
-        surfaceScalarField alphaPhi2BD
-        (
-            upwind<scalar>(mesh, phi).flux(alpha2)
-        );
-
-        // Calculate the flux correction for alpha2
-        alphaPhi2 -= alphaPhi2BD;
-
-        // Further limit the limiter for alpha2
-        if (LTS)
-        {
-            const volScalarField& rDeltaT =
-                fv::localEulerDdt::localRDeltaT(mesh);
-
-            MULES::limiter
-            (
-                lambda,
-                rDeltaT,
-                geometricOneField(),
-                alpha2,
-                alphaPhi2BD,
-                alphaPhi2,
-                zeroField(),
-                zeroField(),
-                oneField(),
-                zeroField()
-            );
-        }
-        else
-        {
-            MULES::limiter
-            (
-                lambda,
-                1.0/runTime.deltaT().value(),
-                geometricOneField(),
-                alpha2,
-                alphaPhi2BD,
-                alphaPhi2,
-                zeroField(),
-                zeroField(),
-                oneField(),
-                zeroField()
-            );
-        }
-
-        // Construct the limited fluxes
-        alphaPhi1 = alphaPhi1BD + lambda*alphaPhi1;
-        alphaPhi2 = alphaPhi2BD + lambda*alphaPhi2;
-
-        // Solve for alpha1
-        solve(fvm::ddt(alpha1) + fvc::div(alphaPhi1));
-
-        // Create the diffusion coefficients for alpha2<->alpha3
-        volScalarField Dc23(D23*max(alpha3, scalar(0))*pos0(alpha2));
-        volScalarField Dc32(D23*max(alpha2, scalar(0))*pos0(alpha3));
-
-        // Add the diffusive flux for alpha3->alpha2
-        alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
-
-        // Solve for alpha2
-        fvScalarMatrix alpha2Eqn
-        (
-            fvm::ddt(alpha2)
-          + fvc::div(alphaPhi2)
-          - fvm::laplacian(Dc23 + Dc32, alpha2)
-        );
-        alpha2Eqn.solve();
-
-        // Construct the complete mass flux
-        rhoPhi =
-              alphaPhi1*(rho1 - rho3)
-            + (alphaPhi2 + alpha2Eqn.flux())*(rho2 - rho3)
-            + phi*rho3;
-
-        alpha3 = 1.0 - alpha1 - alpha2;
-    }
-
-    Info<< "Air phase volume fraction = "
-        << alpha1.weightedAverage(mesh.V()).value()
-        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
-        << endl;
-
-    Info<< "Liquid phase volume fraction = "
-        << alpha2.weightedAverage(mesh.V()).value()
-        << "  Min(" << alpha2.name() << ") = " << min(alpha2).value()
-        << "  Max(" << alpha2.name() << ") = " << max(alpha2).value()
-        << endl;
-}

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnSubCycle.H

@@ -1,43 +0,0 @@
-#include "alphaControls.H"
-
-if (nAlphaSubCycles > 1)
-{
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-    surfaceScalarField rhoPhiSum
-    (
-        IOobject
-        (
-            "rhoPhiSum",
-            runTime.name(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(rhoPhi.dimensions(), 0)
-    );
-
-    for
-    (
-        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-        !(++alphaSubCycle).end();
-    )
-    {
-        #include "alphaEqn.H"
-        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
-    }
-
-    rhoPhi = rhoPhiSum;
-}
-else
-{
-    #include "alphaEqn.H"
-}
-
-{
-    volScalarField rhoNew(alpha1*rho1 + alpha2*rho2 + alpha3*rho3);
-
-    // solve(fvm::ddt(rho) + fvc::div(rhoPhi));
-    // Info<< "density error = "
-    //    << max((mag(rho - rhoNew)/mag(rhoNew))().primitiveField()) << endl;
-
-    rho == rhoNew;
-}

applications/solvers/multiphase/interFoam/interMixingFoam/correctPhi.H

@@ -1,21 +0,0 @@
-// Calculate absolute flux
-// from the mapped surface velocity
-phi = mesh.Sf() & Uf();
-
-correctUphiBCs(U, phi, true);
-
-CorrectPhi
-(
-    phi,
-    U,
-    p_rgh,
-    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-    geometricZeroField(),
-    pressureReference,
-    pimple
-);
-
-#include "continuityErrs.H"
-
-// Make the flux relative to the mesh motion
-fvc::makeRelative(phi, U);

applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H

@@ -1,121 +0,0 @@
-#include "createRDeltaT.H"
-
-Info<< "Reading field p_rgh\n" << endl;
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "createPhi.H"
-
-immiscibleIncompressibleThreePhaseMixture mixture(U);
-
-volScalarField& alpha1(mixture.alpha1());
-volScalarField& alpha2(mixture.alpha2());
-volScalarField& alpha3(mixture.alpha3());
-
-const dimensionedScalar& rho1 = mixture.rho1();
-const dimensionedScalar& rho2 = mixture.rho2();
-const dimensionedScalar& rho3 = mixture.rho3();
-
-dimensionedScalar D23("D23", dimViscosity, mixture);
-
-// Need to store rho for ddt(rho, U)
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.name(),
-        mesh,
-        IOobject::READ_IF_PRESENT
-    ),
-    alpha1*rho1 + alpha2*rho2 + alpha3*rho3
-);
-rho.oldTime();
-
-
-// Mass flux
-// Initialisation does not matter because rhoPhi is reset after the
-// alpha solution before it is used in the U equation.
-surfaceScalarField rhoPhi
-(
-    IOobject
-    (
-        "rhoPhi",
-        runTime.name(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    rho1*phi
-);
-
-
-// Construct incompressible turbulence model
-autoPtr<incompressible::momentumTransportModel> turbulence
-(
-    incompressible::momentumTransportModel::New(U, phi, mixture)
-);
-
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-volScalarField p
-(
-    IOobject
-    (
-        "p",
-        runTime.name(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    p_rgh + rho*gh
-);
-
-pressureReference pressureReference(p, p_rgh, pimple.dict());
-
-if (p_rgh.needReference())
-{
-    p += dimensionedScalar
-    (
-        "p",
-        p.dimensions(),
-        pressureReference.refValue()
-      - getRefCellValue(p, pressureReference.refCell())
-    );
-    p_rgh = p - rho*gh;
-}
-
-mesh.schemes().setFluxRequired(p_rgh.name());
-mesh.schemes().setFluxRequired(alpha2.name());
-
-#include "createMRF.H"
-#include "createFvModels.H"
-#include "createFvConstraints.H"

applications/solvers/multiphase/interFoam/interMixingFoam/immiscibleIncompressibleThreePhaseMixture/immiscibleIncompressibleThreePhaseMixture.H

@@ -1,87 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 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/>.
-
-Class
-    Foam::immiscibleIncompressibleThreePhaseMixture
-
-Description
-    An immiscible incompressible two-phase mixture transport model
-
-SourceFiles
-    immiscibleIncompressibleThreePhaseMixture.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef immiscibleIncompressibleThreePhaseMixture_H
-#define immiscibleIncompressibleThreePhaseMixture_H
-
-#include "incompressibleThreePhaseMixture.H"
-#include "threePhaseInterfaceProperties.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-          Class immiscibleIncompressibleThreePhaseMixture Declaration
-\*---------------------------------------------------------------------------*/
-
-class immiscibleIncompressibleThreePhaseMixture
-:
-    public incompressibleThreePhaseMixture,
-    public threePhaseInterfaceProperties
-{
-public:
-
-    // Constructors
-
-        //- Construct from the mixture velocity
-        immiscibleIncompressibleThreePhaseMixture(const volVectorField& U);
-
-
-    //- Destructor
-    virtual ~immiscibleIncompressibleThreePhaseMixture()
-    {}
-
-
-    // Member Functions
-
-        //- Correct the transport and interface properties
-        virtual void correct()
-        {
-            incompressibleThreePhaseMixture::correct();
-            threePhaseInterfaceProperties::correct();
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.C

@@ -1,200 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "incompressibleThreePhaseMixture.H"
-#include "addToRunTimeSelectionTable.H"
-#include "surfaceFields.H"
-#include "fvc.H"
-
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
-void Foam::incompressibleThreePhaseMixture::calcNu()
-{
-    nuModel1_->correct();
-    nuModel2_->correct();
-    nuModel3_->correct();
-
-    // Average kinematic viscosity calculated from dynamic viscosity
-    nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::incompressibleThreePhaseMixture::incompressibleThreePhaseMixture
-(
-    const fvMesh& mesh
-)
-:
-    IOdictionary
-    (
-        IOobject
-        (
-            "phaseProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    ),
-
-    phase1Name_(wordList(lookup("phases"))[0]),
-    phase2Name_(wordList(lookup("phases"))[1]),
-    phase3Name_(wordList(lookup("phases"))[2]),
-
-    alpha1_
-    (
-        IOobject
-        (
-            IOobject::groupName("alpha", phase1Name_),
-            mesh.time().name(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    ),
-
-    alpha2_
-    (
-        IOobject
-        (
-            IOobject::groupName("alpha", phase2Name_),
-            mesh.time().name(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    ),
-
-    alpha3_
-    (
-        IOobject
-        (
-            IOobject::groupName("alpha", phase3Name_),
-            mesh.time().name(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    ),
-
-    nu_
-    (
-        IOobject
-        (
-            "nu",
-            mesh.time().name(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(dimensionSet(0, 2, -1, 0, 0), 0),
-        calculatedFvPatchScalarField::typeName
-    ),
-
-    nuModel1_(viscosityModel::New(mesh, phase1Name_)),
-    nuModel2_(viscosityModel::New(mesh, phase2Name_)),
-    nuModel3_(viscosityModel::New(mesh, phase3Name_)),
-
-    rho1_("rho", dimDensity, nuModel1_()),
-    rho2_("rho", dimDensity, nuModel2_()),
-    rho3_("rho", dimDensity, nuModel3_())
-{
-    alpha3_ == 1.0 - alpha1_ - alpha2_;
-    calcNu();
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField>
-Foam::incompressibleThreePhaseMixture::mu() const
-{
-    return volScalarField::New
-    (
-        "mu",
-        alpha1_*rho1_*nuModel1_->nu()
-      + alpha2_*rho2_*nuModel2_->nu()
-      + alpha3_*rho3_*nuModel3_->nu()
-    );
-}
-
-
-Foam::tmp<Foam::surfaceScalarField>
-Foam::incompressibleThreePhaseMixture::muf() const
-{
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
-    surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
-    surfaceScalarField alpha3f(fvc::interpolate(alpha3_));
-
-    return surfaceScalarField::New
-    (
-        "mu",
-        alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
-      + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
-      + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
-    );
-}
-
-
-Foam::tmp<Foam::surfaceScalarField>
-Foam::incompressibleThreePhaseMixture::nuf() const
-{
-    surfaceScalarField alpha1f(fvc::interpolate(alpha1_));
-    surfaceScalarField alpha2f(fvc::interpolate(alpha2_));
-    surfaceScalarField alpha3f(fvc::interpolate(alpha3_));
-
-    return surfaceScalarField::New
-    (
-        "nu",
-        (
-            alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
-          + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
-          + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
-        )/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_)
-    );
-}
-
-
-bool Foam::incompressibleThreePhaseMixture::read()
-{
-    if (regIOobject::read())
-    {
-        nuModel1_->lookup("rho") >> rho1_;
-        nuModel2_->lookup("rho") >> rho2_;
-        nuModel3_->lookup("rho") >> rho3_;
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/incompressibleThreePhaseMixture/incompressibleThreePhaseMixture.H

@@ -1,220 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 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/>.
-
-Class
-    Foam::incompressibleThreePhaseMixture
-
-Description
-
-SourceFiles
-    incompressibleThreePhaseMixture.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef incompressibleThreePhaseMixture_H
-#define incompressibleThreePhaseMixture_H
-
-#include "viscosity.H"
-#include "IOdictionary.H"
-#include "viscosityModel.H"
-#include "dimensionedScalar.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-               Class incompressibleThreePhaseMixture Declaration
-\*---------------------------------------------------------------------------*/
-
-class incompressibleThreePhaseMixture
-:
-    public IOdictionary,
-    public viscosity
-{
-    // Private Data
-
-        word phase1Name_;
-        word phase2Name_;
-        word phase3Name_;
-
-        volScalarField alpha1_;
-        volScalarField alpha2_;
-        volScalarField alpha3_;
-
-        volScalarField nu_;
-
-        autoPtr<viscosityModel> nuModel1_;
-        autoPtr<viscosityModel> nuModel2_;
-        autoPtr<viscosityModel> nuModel3_;
-
-        dimensionedScalar rho1_;
-        dimensionedScalar rho2_;
-        dimensionedScalar rho3_;
-
-
-    // Private Member Functions
-
-        //- Calculate and return the laminar viscosity
-        void calcNu();
-
-
-public:
-
-    // Constructors
-
-        //- Construct from a mesh
-        incompressibleThreePhaseMixture(const fvMesh& mesh);
-
-
-    //- Destructor
-    ~incompressibleThreePhaseMixture()
-    {}
-
-
-    // Member Functions
-
-        const word phase1Name() const
-        {
-            return phase1Name_;
-        }
-
-        const word phase2Name() const
-        {
-            return phase2Name_;
-        }
-
-        const word phase3Name() const
-        {
-            return phase3Name_;
-        }
-
-        const volScalarField& alpha1() const
-        {
-            return alpha1_;
-        }
-
-        volScalarField& alpha1()
-        {
-            return alpha1_;
-        }
-
-        const volScalarField& alpha2() const
-        {
-            return alpha2_;
-        }
-
-        volScalarField& alpha2()
-        {
-            return alpha2_;
-        }
-
-        const volScalarField& alpha3() const
-        {
-            return alpha3_;
-        }
-
-        volScalarField& alpha3()
-        {
-            return alpha3_;
-        }
-
-        //- Return const-access to phase1 density
-        const dimensionedScalar& rho1() const
-        {
-            return rho1_;
-        }
-
-        //- Return const-access to phase2 density
-        const dimensionedScalar& rho2() const
-        {
-            return rho2_;
-        };
-
-        //- Return const-access to phase3 density
-        const dimensionedScalar& rho3() const
-        {
-            return rho3_;
-        };
-
-        //- Return const-access to phase1 viscosityModel
-        const viscosityModel& nuModel1() const
-        {
-            return nuModel1_();
-        }
-
-        //- Return const-access to phase2 viscosityModel
-        const viscosityModel& nuModel2() const
-        {
-            return nuModel2_();
-        }
-
-        //- Return const-access to phase3 viscosityModel
-        const viscosityModel& nuModel3() const
-        {
-            return nuModel3_();
-        }
-
-        //- Return the dynamic laminar viscosity
-        tmp<volScalarField> mu() const;
-
-        //- Return the face-interpolated dynamic laminar viscosity
-        tmp<surfaceScalarField> muf() const;
-
-        //- Return the kinematic laminar viscosity
-        tmp<volScalarField> nu() const
-        {
-            return nu_;
-        }
-
-        //- Return the laminar viscosity for patch
-        tmp<scalarField> nu(const label patchi) const
-        {
-            return nu_.boundaryField()[patchi];
-        }
-
-        //- Return the face-interpolated dynamic laminar viscosity
-        tmp<surfaceScalarField> nuf() const;
-
-        //- Correct the laminar viscosity
-        void correct()
-        {
-            calcNu();
-        }
-
-        //- Read base phaseProperties dictionary
-        bool read();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/initCorrectPhi.H

@@ -1,50 +0,0 @@
-tmp<volScalarField> rAU;
-
-if (correctPhi)
-{
-    rAU = new volScalarField
-    (
-        IOobject
-        (
-            "rAU",
-            runTime.name(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar(dimTime/dimDensity, 1)
-    );
-
-    correctUphiBCs(U, phi, true);
-
-    CorrectPhi
-    (
-        phi,
-        U,
-        p_rgh,
-        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-        geometricZeroField(),
-        pressureReference,
-        pimple
-    );
-
-    #include "continuityErrs.H"
-}
-else
-{
-    correctUphiBCs(U, phi, true);
-
-    CorrectPhi
-    (
-        phi,
-        U,
-        p_rgh,
-        dimensionedScalar(dimTime/rho.dimensions(), 1),
-        geometricZeroField(),
-        pressureReference,
-        pimple
-    );
-
-    #include "continuityErrs.H"
-}

applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C

@@ -1,169 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 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
-    interMixingFoam
-
-Description
-    Solver for 3 incompressible fluids, two of which are miscible, using a VOF
-    method to capture the interface, with optional mesh motion and mesh topology
-    changes including adaptive re-meshing.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "interfaceCompression.H"
-#include "CMULES.H"
-#include "localEulerDdtScheme.H"
-#include "subCycle.H"
-#include "immiscibleIncompressibleThreePhaseMixture.H"
-#include "incompressibleMomentumTransportModels.H"
-#include "pimpleControl.H"
-#include "pressureReference.H"
-#include "fvModels.H"
-#include "fvConstraints.H"
-#include "CorrectPhi.H"
-#include "fvcSmooth.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "initContinuityErrs.H"
-    #include "createDyMControls.H"
-    #include "createFields.H"
-    #include "initCorrectPhi.H"
-    #include "createUfIfPresent.H"
-
-    turbulence->validate();
-
-    if (!LTS)
-    {
-        #include "readTimeControls.H"
-        #include "CourantNo.H"
-        #include "setInitialDeltaT.H"
-    }
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (pimple.run(runTime))
-    {
-        #include "readDyMControls.H"
-
-        if (LTS)
-        {
-            #include "setRDeltaT.H"
-        }
-        else
-        {
-            #include "CourantNo.H"
-            #include "alphaCourantNo.H"
-            #include "setDeltaT.H"
-        }
-
-        fvModels.preUpdateMesh();
-
-        // Update the mesh for topology change, mesh to mesh mapping
-        mesh.update();
-
-        runTime++;
-
-        Info<< "Time = " << runTime.userTimeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            if (pimple.firstPimpleIter() || pimple.moveMeshOuterCorrectors())
-            {
-                // Move the mesh
-                mesh.move();
-
-                if (mesh.changing())
-                {
-                    gh = (g & mesh.C()) - ghRef;
-                    ghf = (g & mesh.Cf()) - ghRef;
-
-                    MRF.update();
-
-                    if (correctPhi)
-                    {
-                        #include "correctPhi.H"
-                    }
-
-                    mixture.correct();
-
-                    if (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
-                }
-            }
-
-            fvModels.correct();
-
-            #include "alphaEqnSubCycle.H"
-
-            mixture.correct();
-
-            if (pimple.predictTransport())
-            {
-                turbulence->predict();
-            }
-
-            #include "UEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.correctTransport())
-            {
-                turbulence->correct();
-            }
-        }
-
-        #include "continuityErrs.H"
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/pEqn.H

@@ -1,94 +0,0 @@
-{
-    if (rAU.valid())
-    {
-        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
-        (
-            pressureReference.refCell(),
-            getRefCellValue(p_rgh, pressureReference.refCell())
-        );
-
-        p_rghEqn.solve();
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA - p_rghEqn.flux();
-
-            p_rgh.relax();
-
-            U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
-            U.correctBoundaryConditions();
-            fvConstraints.constrain(U);
-        }
-    }
-
-    #include "continuityErrs.H"
-
-    // Correct Uf if the mesh is moving
-    fvc::correctUf(Uf, U, phi, MRF);
-
-    // 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(),
-            pressureReference.refValue()
-          - getRefCellValue(p, pressureReference.refCell())
-        );
-        p_rgh = p - rho*gh;
-    }
-
-    if (!correctPhi)
-    {
-        rAU.clear();
-    }
-}

applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C

@@ -1,218 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "threePhaseInterfaceProperties.H"
-#include "alphaContactAngleFvPatchScalarField.H"
-#include "unitConversion.H"
-#include "surfaceInterpolate.H"
-#include "fvcDiv.H"
-#include "fvcGrad.H"
-#include "fvcSnGrad.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::threePhaseInterfaceProperties::correctContactAngle
-(
-    surfaceVectorField::Boundary& nHatb
-) const
-{
-    const volScalarField::Boundary& alpha1Bf = alpha1_.boundaryField();
-    const volScalarField::Boundary& alpha2Bf = alpha2_.boundaryField();
-    const volScalarField::Boundary& alpha3Bf = alpha3_.boundaryField();
-    const volVectorField::Boundary& UBf = U_.boundaryField();
-
-    const fvMesh& mesh = alpha1_.mesh();
-    const fvBoundaryMesh& boundary = mesh.boundary();
-
-    forAll(boundary, patchi)
-    {
-        if (isA<alphaContactAngleFvPatchScalarField>(alpha1Bf[patchi]))
-        {
-            const alphaContactAngleFvPatchScalarField& a2cap =
-                refCast<const alphaContactAngleFvPatchScalarField>
-                (alpha2Bf[patchi]);
-
-            const alphaContactAngleFvPatchScalarField& a3cap =
-                refCast<const alphaContactAngleFvPatchScalarField>
-                (alpha3Bf[patchi]);
-
-            scalarField twoPhaseAlpha2(max(a2cap, scalar(0)));
-            scalarField twoPhaseAlpha3(max(a3cap, scalar(0)));
-
-            scalarField sumTwoPhaseAlpha
-            (
-                twoPhaseAlpha2 + twoPhaseAlpha3 + small
-            );
-
-            twoPhaseAlpha2 /= sumTwoPhaseAlpha;
-            twoPhaseAlpha3 /= sumTwoPhaseAlpha;
-
-            fvsPatchVectorField& nHatp = nHatb[patchi];
-
-            scalarField theta
-            (
-                degToRad
-                (
-                   twoPhaseAlpha2*(180 - a2cap.theta(UBf[patchi], nHatp))
-                 + twoPhaseAlpha3*(180 - a3cap.theta(UBf[patchi], nHatp))
-                )
-            );
-
-            vectorField nf(boundary[patchi].nf());
-
-            // Reset nHatPatch to correspond to the contact angle
-
-            scalarField a12(nHatp & nf);
-
-            scalarField b1(cos(theta));
-
-            scalarField b2(nHatp.size());
-
-            forAll(b2, facei)
-            {
-                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
-            }
-
-            scalarField det(1.0 - a12*a12);
-
-            scalarField a((b1 - a12*b2)/det);
-            scalarField b((b2 - a12*b1)/det);
-
-            nHatp = a*nf + b*nHatp;
-
-            nHatp /= (mag(nHatp) + deltaN_.value());
-        }
-    }
-}
-
-
-void Foam::threePhaseInterfaceProperties::calculateK()
-{
-    const fvMesh& mesh = alpha1_.mesh();
-    const surfaceVectorField& Sf = mesh.Sf();
-
-    // Cell gradient of alpha
-    volVectorField gradAlpha(fvc::grad(alpha1_));
-
-    // Interpolated face-gradient of alpha
-    surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
-
-    // Face unit interface normal
-    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
-
-    correctContactAngle(nHatfv.boundaryFieldRef());
-
-    // Face unit interface normal flux
-    nHatf_ = nHatfv & Sf;
-
-    // Simple expression for curvature
-    K_ = -fvc::div(nHatf_);
-
-    // Complex expression for curvature.
-    // Correction is formally zero but numerically non-zero.
-    // volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
-    // nHat.boundaryField() = nHatfv.boundaryField();
-    // K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
-(
-    const IOdictionary& dict,
-    volScalarField& alpha1,
-    volScalarField& alpha2,
-    volScalarField& alpha3,
-    const volVectorField& U
-)
-:
-    alpha1_(alpha1),
-    alpha2_(alpha2),
-    alpha3_(alpha3),
-    U_(U),
-
-    cAlpha_
-    (
-        U.mesh().solution().solverDict(alpha1_.name()).lookup<scalar>("cAlpha")
-    ),
-    sigma12_("sigma12", dimensionSet(1, 0, -2, 0, 0), dict),
-    sigma13_("sigma13", dimensionSet(1, 0, -2, 0, 0), dict),
-
-    deltaN_
-    (
-        "deltaN",
-        1e-8/pow(average(U.mesh().V()), 1.0/3.0)
-    ),
-
-    nHatf_
-    (
-        IOobject
-        (
-            "nHatf",
-            alpha1_.time().name(),
-            alpha1_.mesh()
-        ),
-        alpha1_.mesh(),
-        dimensionedScalar(dimArea, 0)
-    ),
-
-    K_
-    (
-        IOobject
-        (
-            "interfaceProperties:K",
-            alpha1_.time().name(),
-            alpha1_.mesh()
-        ),
-        alpha1_.mesh(),
-        dimensionedScalar(dimless/dimLength, 0)
-    )
-{
-    calculateK();
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::surfaceScalarField>
-Foam::threePhaseInterfaceProperties::surfaceTensionForce() const
-{
-    return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::threePhaseInterfaceProperties::nearInterface() const
-{
-    return max
-    (
-        pos0(alpha1_ - 0.01)*pos0(0.99 - alpha1_),
-        pos0(alpha2_ - 0.01)*pos0(0.99 - alpha2_)
-    );
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H

@@ -1,192 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 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/>.
-
-Class
-    Foam::threePhaseInterfaceProperties
-
-Description
-    Properties to aid interFoam :
-    1. Correct the alpha boundary condition for dynamic contact angle.
-    2. Calculate interface curvature.
-
-SourceFiles
-    threePhaseInterfaceProperties.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef threePhaseInterfaceProperties_H
-#define threePhaseInterfaceProperties_H
-
-#include "incompressibleThreePhaseMixture.H"
-#include "surfaceFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                Class threePhaseInterfaceProperties Declaration
-\*---------------------------------------------------------------------------*/
-
-class threePhaseInterfaceProperties
-{
-    // Private Data
-
-        //- Volume fraction of phase 1
-        volScalarField& alpha1_;
-
-        //- Volume fraction of phase 2
-        volScalarField& alpha2_;
-
-        //- Volume fraction of phase 3
-        volScalarField& alpha3_;
-
-        //- Mixture velocity
-        const volVectorField& U_;
-
-        //- Compression coefficient
-        scalar cAlpha_;
-
-        //- Surface tension 1-2
-        dimensionedScalar sigma12_;
-
-        //- Surface tension 1-3
-        dimensionedScalar sigma13_;
-
-        //- Stabilisation for normalisation of the interface normal
-        const dimensionedScalar deltaN_;
-
-        //- ...
-        surfaceScalarField nHatf_;
-
-        //- ...
-        volScalarField K_;
-
-
-    // Private Member Functions
-
-        //- Correction for the boundary condition on the unit normal nHat on
-        //  walls to produce the correct contact dynamic angle.
-        //  Calculated from the component of U parallel to the wall
-        void correctContactAngle
-        (
-            surfaceVectorField::Boundary& nHat
-        ) const;
-
-        //- Re-calculate the interface curvature
-        void calculateK();
-
-
-public:
-
-    // Constructors
-
-        //- Construct from dictionary, volume fraction fields and mixture
-        //  velocity
-        threePhaseInterfaceProperties
-        (
-            const IOdictionary& dict,
-            volScalarField& alpha1,
-            volScalarField& alpha2,
-            volScalarField& alpha3,
-            const volVectorField& U
-        );
-
-        //- Disallow default bitwise copy construction
-        threePhaseInterfaceProperties
-        (
-            const threePhaseInterfaceProperties&
-        ) = delete;
-
-
-    // Member Functions
-
-        inline const volVectorField& U() const
-        {
-            return U_;
-        }
-
-        scalar cAlpha() const
-        {
-            return cAlpha_;
-        }
-
-        const dimensionedScalar& deltaN() const
-        {
-            return deltaN_;
-        }
-
-        const surfaceScalarField& nHatf() const
-        {
-            return nHatf_;
-        }
-
-        const volScalarField& K() const
-        {
-            return K_;
-        }
-
-        tmp<volScalarField> sigma() const
-        {
-            volScalarField limitedAlpha2(max(alpha2_, scalar(0)));
-            volScalarField limitedAlpha3(max(alpha3_, scalar(0)));
-
-            return
-                (limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_)
-               /(limitedAlpha2 + limitedAlpha3 + small);
-        }
-
-        tmp<volScalarField> sigmaK() const
-        {
-            return sigma()*K_;
-        }
-
-        tmp<surfaceScalarField> surfaceTensionForce() const;
-
-        //- Indicator of the proximity of the interface
-        //  Field values are 1 near and 0 away for the interface.
-        tmp<volScalarField> nearInterface() const;
-
-        void correct()
-        {
-            calculateK();
-        }
-
-
-    // Member Operators
-
-        //- Disallow default bitwise assignment
-        void operator=(const threePhaseInterfaceProperties&) = delete;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/multiphase/interFoam/pEqn.H

@@ -1,110 +0,0 @@
-{
-    if (rAU.valid())
-    {
-        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);
-
-    // Cache the phase change pressure source
-    fvScalarMatrix Sp_rgh
-    (
-        fvModels.source
-        (
-            volScalarField::New
-            (
-                "1",
-                mesh,
-                dimensionedScalar(dimless/dimPressure, 1)
-            ),
-            p_rgh
-        )
-    );
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix p_rghEqn
-        (
-            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-         == Sp_rgh
-        );
-
-        p_rghEqn.setReference
-        (
-            pressureReference.refCell(),
-            getRefCellValue(p_rgh, pressureReference.refCell())
-        );
-
-        p_rghEqn.solve();
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + p_rghEqn.flux();
-
-            p_rgh.relax();
-
-            U = HbyA + rAU()*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
-            U.correctBoundaryConditions();
-            fvConstraints.constrain(U);
-        }
-    }
-
-    #include "continuityErrs.H"
-
-    // Correct Uf if the mesh is moving
-    fvc::correctUf(Uf, U, phi, MRF);
-
-    // 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(),
-            pressureReference.refValue()
-          - getRefCellValue(p, pressureReference.refCell())
-        );
-        p_rgh = p - rho*gh;
-    }
-
-    if (!correctPhi)
-    {
-        rAU.clear();
-    }
-}

applications/solvers/multiphase/interFoam/rhofs.H

@@ -1,2 +0,0 @@
-const dimensionedScalar& rho1f(rho1);
-const dimensionedScalar& rho2f(rho2);

bin/interFoam

@@ -0,0 +1,46 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     | Website:  https://openfoam.org
+#   \\  /    A nd           | Copyright (C) 2022 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/>.
+#
+# Script
+#     interFoam
+#
+# Description
+#     Script to inform the user that interFoam has been superseded and replaced
+#     by the more general incompressibleVoF solver module executed by the
+#     foamRun application.
+#
+#------------------------------------------------------------------------------
+
+cat <<EOF
+
+interFoam has been superseded and replaced by the more general
+incompressibleVoF solver module executed by the foamRun application:
+
+    foamRun -solver incompressibleVoF
+
+EOF
+
+exec env foamRun -solver incompressibleVoF "$@"
+
+#------------------------------------------------------------------------------

etc/config.sh/bash_completion

@@ -1819,56 +1819,6 @@ _insideCells_ ()
 }
 complete -o filenames -o nospace -F _insideCells_ insideCells
 
-_interFoam_ ()
-{
-    local cur="${COMP_WORDS[COMP_CWORD]}"
-    local prev="${COMP_WORDS[COMP_CWORD-1]}"
-    local line=${COMP_LINE}
-    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
-
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvConstraints -listFvModels -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
-    for o in $used ; do opts="${opts/$o/}" ; done
-    extra=""
-
-    [ "$COMP_CWORD" = 1 ] || \
-    case "$prev" in
-        -case)
-            opts="" ; extra="-d" ;;
-        -fileHandler)
-            opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
-            opts="" ; extra="" ;;
-       *) ;;
-    esac
-    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
-}
-complete -o filenames -o nospace -F _interFoam_ interFoam
-
-_interMixingFoam_ ()
-{
-    local cur="${COMP_WORDS[COMP_CWORD]}"
-    local prev="${COMP_WORDS[COMP_CWORD-1]}"
-    local line=${COMP_LINE}
-    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
-
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvConstraints -listFvModels -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
-    for o in $used ; do opts="${opts/$o/}" ; done
-    extra=""
-
-    [ "$COMP_CWORD" = 1 ] || \
-    case "$prev" in
-        -case)
-            opts="" ; extra="-d" ;;
-        -fileHandler)
-            opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
-            opts="" ; extra="" ;;
-       *) ;;
-    esac
-    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
-}
-complete -o filenames -o nospace -F _interMixingFoam_ interMixingFoam
-
 _kivaToFoam_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"
@@ -5086,6 +5036,31 @@ _foamTags_ ()
 }
 complete -o filenames -o nospace -F _foamTags_ foamTags
 
+_interFoam_ ()
+{
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    local line=${COMP_LINE}
+    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
+
+    opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -parallel -roots -solver -srcDoc"
+    for o in $used ; do opts="${opts/$o/}" ; done
+    extra=""
+
+    [ "$COMP_CWORD" = 1 ] || \
+    case "$prev" in
+        -case)
+            opts="" ; extra="-d" ;;
+        -fileHandler)
+            opts="uncollated collated masterUncollated" ; extra="" ;;
+        -hostRoots|-libs|-roots|-solver)
+            opts="" ; extra="" ;;
+       *) ;;
+    esac
+    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
+}
+complete -o filenames -o nospace -F _interFoam_ interFoam
+
 _mpirunDebug_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"

src/twoPhaseModels/Allwmake

@@ -13,8 +13,6 @@ wmake $targetType compressibleInterfaceProperties
 wmake $targetType twoPhaseMixture
 wmake $targetType incompressibleTwoPhaseMixture
 wmake $targetType compressibleTwoPhaseMixture
-wmake $targetType immiscibleIncompressibleTwoPhaseMixture
-wmake $targetType immiscibleCompressibleTwoPhaseMixture
 
 wmake $targetType incompressibleCavitation
 wmake $targetType compressibleCavitation

src/twoPhaseModels/VoF/alphaScheme.H

@@ -1,20 +0,0 @@
-static const word divAlphaName("div(phi,alpha)");
-
-const word alphaScheme(mesh.schemes().div(divAlphaName)[1].wordToken());
-
-ITstream compressionScheme
-(
-    compressionSchemes.found(alphaScheme)
-  ? mesh.schemes().div(divAlphaName)
-  : ITstream
-    (
-        divAlphaName,
-        tokenList
-        {
-            word("Gauss"),
-            word("interfaceCompression"),
-            alphaScheme,
-            alphaControls.lookup<scalar>("cAlpha")
-        }
-    )
-);

src/twoPhaseModels/VoF/createAlphaFluxes.H

@@ -1,25 +0,0 @@
-typeIOobject<surfaceScalarField> alphaPhi1Header
-(
-    IOobject::groupName("alphaPhi", alpha1.group()),
-    runTime.name(),
-    mesh,
-    IOobject::READ_IF_PRESENT,
-    IOobject::AUTO_WRITE
-);
-
-const bool alphaRestart = alphaPhi1Header.headerOk();
-
-if (alphaRestart)
-{
-    Info << "Restarting alpha" << endl;
-}
-
-// Phase-1 volumetric flux
-surfaceScalarField alphaPhi1
-(
-    alphaPhi1Header,
-    phi*fvc::interpolate(alpha1)
-);
-
-// MULES Correction
-tmp<surfaceScalarField> talphaPhi1Corr0;

src/twoPhaseModels/VoF/setRDeltaT.H

@@ -1,147 +0,0 @@
-{
-    volScalarField& rDeltaT = trDeltaT.ref();
-
-    const dictionary& pimpleDict = pimple.dict();
-
-    const scalar maxCo
-    (
-        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
-    );
-
-    const scalar maxAlphaCo
-    (
-        pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
-    );
-
-    const scalar rDeltaTSmoothingCoeff
-    (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
-    );
-
-    const label nAlphaSpreadIter
-    (
-        pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
-    );
-
-    const scalar alphaSpreadDiff
-    (
-        pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
-    );
-
-    const scalar alphaSpreadMax
-    (
-        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
-    );
-
-    const scalar alphaSpreadMin
-    (
-        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
-    );
-
-    const label nAlphaSweepIter
-    (
-        pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
-    );
-
-    const scalar maxDeltaT
-    (
-        pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great)
-    );
-
-    const scalar minDeltaT
-    (
-        pimpleDict.lookupOrDefault<scalar>("minDeltaT", small)
-    );
-
-    const volScalarField rDeltaT0("rDeltaT0", rDeltaT);
-
-    // Set the reciprocal time-step from the local Courant number
-    // and maximum and minimum time-steps
-    rDeltaT.ref() = min
-    (
-        1/dimensionedScalar(dimTime, minDeltaT),
-        max
-        (
-            1/dimensionedScalar(dimTime, maxDeltaT),
-            fvc::surfaceSum(mag(phi))()()
-           /((2*maxCo)*mesh.V())
-        )
-    );
-
-    if (maxAlphaCo < maxCo)
-    {
-        // Further limit the reciprocal time-step
-        // in the vicinity of the interface
-
-        volScalarField alpha1Bar(fvc::average(alpha1));
-
-        rDeltaT.ref() = max
-        (
-            rDeltaT(),
-            pos0(alpha1Bar() - alphaSpreadMin)
-           *pos0(alphaSpreadMax - alpha1Bar())
-           *fvc::surfaceSum(mag(phi))()()
-           /((2*maxAlphaCo)*mesh.V())
-        );
-    }
-
-    // Update the boundary values of the reciprocal time-step
-    rDeltaT.correctBoundaryConditions();
-
-    Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
-
-    if (rDeltaTSmoothingCoeff < 1.0)
-    {
-        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
-    }
-
-    if (nAlphaSpreadIter > 0)
-    {
-        fvc::spread
-        (
-            rDeltaT,
-            alpha1,
-            nAlphaSpreadIter,
-            alphaSpreadDiff,
-            alphaSpreadMax,
-            alphaSpreadMin
-        );
-    }
-
-    if (nAlphaSweepIter > 0)
-    {
-        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
-    }
-
-    Info<< "Smoothed flow time scale min/max = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
-
-    // Limit rate of change of time scale
-    // - reduce as much as required
-    // - only increase at a fraction of old time scale
-    if
-    (
-        pimpleDict.found("rDeltaTDampingCoeff")
-     && runTime.timeIndex() > runTime.startTimeIndex() + 1
-    )
-    {
-        // Damping coefficient (1-0)
-        const scalar rDeltaTDampingCoeff
-        (
-            pimpleDict.lookup<scalar>("rDeltaTDampingCoeff")
-        );
-
-        rDeltaT = max
-        (
-            rDeltaT,
-            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
-        );
-
-        Info<< "Damped flow time scale min/max = "
-            << gMin(1/rDeltaT.primitiveField())
-            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
-    }
-}

src/twoPhaseModels/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.H

@@ -162,7 +162,7 @@ public:
         }
 
         //- Return mixture density [kg/m^3]
-        const volScalarField& rho() const
+        virtual const volScalarField& rho() const
         {
             return rho_;
         }

src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/Make/files

@@ -1,3 +0,0 @@
-immiscibleCompressibleTwoPhaseMixture.C
-
-LIB = $(FOAM_LIBBIN)/libimmiscibleCompressibleTwoPhaseMixture

src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/Make/options

@@ -1,19 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/physicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/twoPhaseProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-LIB_LIBS = \
-    -lphysicalProperties \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -ltwoPhaseProperties \
-    -linterfaceProperties \
-    -ltwoPhaseMixture \
-    -lcompressibleTwoPhaseMixture \
-    -lfiniteVolume

src/twoPhaseModels/immiscibleCompressibleTwoPhaseMixture/immiscibleCompressibleTwoPhaseMixture.H

@@ -1,87 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 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/>.
-
-Class
-    Foam::immiscibleCompressibleTwoPhaseMixture
-
-Description
-    Class to represent a mixture of two immiscible rhoThermo-based phases
-
-SourceFiles
-    immiscibleCompressibleTwoPhaseMixture.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef immiscibleCompressibleTwoPhaseMixture_H
-#define immiscibleCompressibleTwoPhaseMixture_H
-
-#include "compressibleTwoPhaseMixture.H"
-#include "interfaceProperties.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-           Class immiscibleCompressibleTwoPhaseMixture Declaration
-\*---------------------------------------------------------------------------*/
-
-class immiscibleCompressibleTwoPhaseMixture
-:
-    public compressibleTwoPhaseMixture,
-    public interfaceProperties
-{
-
-public:
-
-    // Constructors
-
-        //- Construct from the mixture velocity
-        immiscibleCompressibleTwoPhaseMixture(const volVectorField& U);
-
-
-    //- Destructor
-    virtual ~immiscibleCompressibleTwoPhaseMixture()
-    {}
-
-
-    // Member Functions
-
-        //- Correct the mixture and interface properties
-        virtual void correct();
-
-        //- Read base phaseProperties dictionary
-        virtual bool read();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/Make/files

@@ -1,3 +0,0 @@
-immiscibleIncompressibleTwoPhaseMixture.C
-
-LIB = $(FOAM_LIBBIN)/libimmiscibleIncompressibleTwoPhaseMixture

src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/Make/options

@@ -1,15 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/physicalProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/twoPhaseProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-LIB_LIBS = \
-    -lphysicalProperties \
-    -ltwoPhaseProperties \
-    -linterfaceProperties \
-    -ltwoPhaseMixture \
-    -lincompressibleTwoPhaseMixture \
-    -lfiniteVolume

src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C

@@ -1,54 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 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/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "immiscibleIncompressibleTwoPhaseMixture.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::immiscibleIncompressibleTwoPhaseMixture::
-immiscibleIncompressibleTwoPhaseMixture(const volVectorField& U)
-:
-    incompressibleTwoPhaseMixture(U.mesh()),
-    interfaceProperties(*this, alpha1(), alpha2(), U)
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-void Foam::immiscibleIncompressibleTwoPhaseMixture::correct()
-{
-    incompressibleTwoPhaseMixture::correct();
-}
-
-
-bool Foam::immiscibleIncompressibleTwoPhaseMixture::read()
-{
-    return
-        incompressibleTwoPhaseMixture::read()
-     && interfaceProperties::read();
-}
-
-
-// ************************************************************************* //