solvers::incompressibleDriftFlux: New solver module for two-phase flow with drift-flux

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

Class
    Foam::solvers::incompressibleDriftFlux

Description
    Solver module for 2 incompressible fluids using the mixture approach with
    the drift-flux approximation for relative motion of the phases, 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 with mixture transport modelling in which a
    single laminar, RAS or LES model is selected to model the momentum stress.

    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
    incompressibleDriftFlux.C

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

applications/solvers/modules/Allwmake

@@ -20,6 +20,7 @@ compressibleVoF/Allwmake $targetType $*
 wmake $targetType multiphaseVoFSolver
 wmake $targetType incompressibleMultiphaseVoF
 wmake $targetType compressibleMultiphaseVoF
+wmake $targetType incompressibleDriftFlux
 multiphaseEuler/Allwmake $targetType $*
 wmake $targetType incompressibleDenseParticleFluid
 isothermalFilm/Allwmake $targetType $*

applications/solvers/modules/incompressibleDriftFlux/Allwclean

@@ -1,8 +1,9 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
 
-wclean libso mixtureViscosityModels
-wclean libso relativeVelocityModels
+wclean mixtureViscosityModels
+wclean $targetType relativeVelocityModels
+
 wclean
 
 #------------------------------------------------------------------------------

applications/solvers/modules/incompressibleDriftFlux/Allwmake

@@ -6,6 +6,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 wmake $targetType mixtureViscosityModels
 wmake $targetType relativeVelocityModels
+
 wmake $targetType
 
 #------------------------------------------------------------------------------

applications/solvers/modules/incompressibleDriftFlux/Make/files

@@ -0,0 +1,6 @@
+incompressibleDriftFluxMixture/incompressibleDriftFluxMixture.C
+alphaSuSp.C
+pressureCorrector.C
+incompressibleDriftFlux.C
+
+LIB = $(FOAM_LIBBIN)/libincompressibleDriftFlux

applications/solvers/modules/incompressibleDriftFlux/Make/options

@@ -0,0 +1,32 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
+    -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)/twoPhaseModels/incompressibleTwoPhases/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lincompressibleDriftFluxViscosityModel \
+    -lincompressibleDriftFluxRelativeVelocityModels \
+    -ltwoPhaseVoFSolver \
+    -lphysicalProperties \
+    -linterfaceCompression \
+    -lincompressibleTwoPhases \
+    -ltwoPhaseProperties \
+    -linterfaceProperties \
+    -lmomentumTransportModels \
+    -lcompressibleMomentumTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvModels \
+    -lfvConstraints \
+    -lsampling

applications/solvers/modules/incompressibleDriftFlux/alphaSuSp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,33 +21,45 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Global
-    CourantNo
-
-Description
-    Calculates and outputs the mean and maximum Courant Numbers.
-
 \*---------------------------------------------------------------------------*/
 
-scalar CoNum = 0.0;
-scalar meanCoNum = 0.0;
+#include "incompressibleDriftFlux.H"
+#include "fvcFlux.H"
+#include "fvcDiv.H"
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::solvers::incompressibleDriftFlux::alphaPhi
+(
+    const surfaceScalarField& phi,
+    const volScalarField& alpha,
+    const dictionary& alphaControls
+)
 {
-    const scalarField sumPhi
+    return fvc::flux
     (
-        fvc::surfaceSum
-        (
-            mag(phi) + mag(fvc::flux(relativeVelocity->Udm()))
-        )().primitiveField()
+        phi + fvc::flux(relativeVelocity->Udm()),
+        alpha,
+        divAlphaName
     );
-
-    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
-
-    meanCoNum =
-        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
 }
 
-Info<< "Courant Number mean: " << meanCoNum
-    << " max: " << CoNum << endl;
+void Foam::solvers::incompressibleDriftFlux::alphaSuSp
+(
+    tmp<volScalarField::Internal>& Su,
+    tmp<volScalarField::Internal>& Sp
+)
+{
+    if (divergent())
+    {
+        // Phase change alpha1 source
+        const fvScalarMatrix alphaSup(fvModels().source(alpha1));
+
+        Su = alphaSup.Su();
+        Sp = alphaSup.Sp();
+    }
+}
+
 
 // ************************************************************************* //

applications/solvers/modules/incompressibleDriftFlux/incompressibleDriftFlux.C

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022-2023 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 "incompressibleDriftFlux.H"
+#include "localEulerDdtScheme.H"
+#include "fvCorrectPhi.H"
+#include "geometricZeroField.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "fvmDdt.H"
+
+#include "fvcDdt.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(incompressibleDriftFlux, 0);
+    addToRunTimeSelectionTable(solver, incompressibleDriftFlux, fvMesh);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::incompressibleDriftFlux::incompressibleDriftFlux(fvMesh& mesh)
+:
+    twoPhaseVoFSolver
+    (
+        mesh,
+        autoPtr<twoPhaseVoFMixture>(new incompressibleDriftFluxMixture(mesh))
+    ),
+
+    mixture
+    (
+        refCast<incompressibleDriftFluxMixture>(twoPhaseVoFSolver::mixture)
+       .initialise(U)
+    ),
+
+    p
+    (
+        IOobject
+        (
+            "p",
+            runTime.name(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p_rgh + rho*buoyancy.gh
+    ),
+
+    pressureReference_
+    (
+        p,
+        p_rgh,
+        pimple.dict()
+    ),
+
+    relativeVelocity
+    (
+        relativeVelocityModel::New(mixture, mixture, buoyancy.g)
+    ),
+
+    momentumTransport
+    (
+        compressible::momentumTransportModel::New(rho, U, rhoPhi, mixture)
+    )
+{
+    // Read the controls
+    readControls();
+
+    if (correctPhi || mesh.topoChanging())
+    {
+        rAU = new volScalarField
+        (
+            IOobject
+            (
+                "rAU",
+                runTime.name(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimTime/dimDensity, 1)
+        );
+    }
+
+    if (!runTime.restart() || !divergent())
+    {
+        correctUphiBCs(U_, phi_, true);
+
+        fv::correctPhi
+        (
+            phi_,
+            U,
+            p_rgh,
+            rAU,
+            autoPtr<volScalarField>(),
+            pressureReference(),
+            pimple
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::incompressibleDriftFlux::~incompressibleDriftFlux()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleDriftFlux::prePredictor()
+{
+    VoFSolver::prePredictor();
+    alphaPredictor();
+
+    // Apply the diffusion term separately to allow implicit solution
+    // and boundedness of the explicit advection
+    {
+        fvScalarMatrix alpha1Eqn
+        (
+            fvm::ddt(alpha1) - fvc::ddt(alpha1)
+          - fvm::laplacian(momentumTransport->nut(), alpha1)
+        );
+
+        alpha1Eqn.solve(alpha1.name() + "Diffusion");
+
+        alphaPhi1 += alpha1Eqn.flux();
+        alpha2 = 1.0 - alpha1;
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+            << endl;
+    }
+
+    mixture.correct();
+
+    const dimensionedScalar& rho1 = mixture.rho1();
+    const dimensionedScalar& rho2 = mixture.rho2();
+
+    // Calculate the mass-flux from the accumulated alphaPhi1
+    rhoPhi = (alphaPhi1*(rho1 - rho2) + phi*rho2);
+
+    relativeVelocity->correct();
+
+    if (pimple.predictTransport())
+    {
+        momentumTransport->predict();
+    }
+}
+
+
+void Foam::solvers::incompressibleDriftFlux::thermophysicalPredictor()
+{}
+
+
+void Foam::solvers::incompressibleDriftFlux::postCorrector()
+{
+    if (pimple.correctTransport())
+    {
+        momentumTransport->correct();
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/incompressibleDriftFlux/incompressibleDriftFlux.H

@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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::incompressibleDriftFlux
+
+Description
+    Solver module for 2 incompressible fluids using the mixture approach with
+    the drift-flux approximation for relative motion of the phases, 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 with mixture transport modelling in which a
+    single laminar, RAS or LES model is selected to model the momentum stress.
+
+    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
+    incompressibleDriftFlux.C
+
+See also
+    Foam::solvers::VoFSolver
+    Foam::solvers::twoPhaseVoFSolver
+    Foam::solvers::compressibleVoF
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef incompressibleDriftFlux_H
+#define incompressibleDriftFlux_H
+
+#include "twoPhaseVoFSolver.H"
+#include "incompressibleDriftFluxMixture.H"
+#include "relativeVelocityModel.H"
+#include "momentumTransportModel.H"
+#include "compressibleMomentumTransportModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class incompressibleDriftFlux Declaration
+\*---------------------------------------------------------------------------*/
+
+class incompressibleDriftFlux
+:
+    public twoPhaseVoFSolver
+{
+
+protected:
+
+    // Phase properties
+
+        //- The compressible two-phase mixture
+        incompressibleDriftFluxMixture& mixture;
+
+
+    // Thermophysical properties
+
+        //- Static pressure field
+        volScalarField p;
+
+
+    // Pressure reference
+
+        //- Pressure reference
+        Foam::pressureReference pressureReference_;
+
+
+    // Momentum transport
+
+        //- Pointer to the dispersed phase relative velocity model
+        autoPtr<relativeVelocityModel> relativeVelocity;
+
+        //- Pointer to the momentum transport model
+        autoPtr<compressible::momentumTransportModel> momentumTransport;
+
+
+    // Protected Member Functions
+
+        //- Return the pressure reference
+        virtual const Foam::pressureReference& pressureReference() const
+        {
+            return pressureReference_;
+        }
+
+        //- The flow is incompressible
+        virtual bool incompressible() const
+        {
+            return true;
+        }
+
+        //- Is the flow divergent?
+        //  i.e. includes phase-fraction sources
+        virtual bool divergent() const
+        {
+            return fvModels().addsSupToField(alpha1.name());
+        }
+
+        //- Return the mixture compressibility/density
+        //  Not required for incompressible fluids
+        virtual tmp<volScalarField> psiByRho() const
+        {
+            return tmp<volScalarField>(nullptr);
+        }
+
+        virtual tmp<surfaceScalarField> alphaPhi
+        (
+            const surfaceScalarField& phi,
+            const volScalarField& alpha,
+            const dictionary& alphaControls
+        );
+
+        //- 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
+                relativeVelocity->divDevTau()
+              + momentumTransport->divDevTau(U);
+        }
+
+
+public:
+
+    //- Runtime type information
+    TypeName("incompressibleDriftFlux");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        incompressibleDriftFlux(fvMesh& mesh);
+
+        //- Disallow default bitwise copy construction
+        incompressibleDriftFlux(const incompressibleDriftFlux&) = delete;
+
+
+    //- Destructor
+    virtual ~incompressibleDriftFlux();
+
+
+    // 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();
+
+        //- Correct the momentum and thermophysical transport modelling
+        virtual void postCorrector();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const incompressibleDriftFlux&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/incompressibleDriftFlux/incompressibleDriftFluxMixture/incompressibleDriftFluxMixture.C

@@ -23,36 +23,31 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "incompressibleTwoPhaseInteractingMixture.H"
-#include "mixtureViscosityModel.H"
+#include "incompressibleDriftFluxMixture.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(incompressibleTwoPhaseInteractingMixture, 0);
+    defineTypeNameAndDebug(incompressibleDriftFluxMixture, 0);
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::incompressibleTwoPhaseInteractingMixture::
-incompressibleTwoPhaseInteractingMixture
+Foam::incompressibleDriftFluxMixture::incompressibleDriftFluxMixture
 (
-    volVectorField& U,
-    const surfaceScalarField& phi
+    const fvMesh& mesh
 )
 :
-    twoPhaseMixture(U.mesh()),
+    twoPhaseVoFMixture(mesh),
 
-    U_(U),
-
-    nucModel_(viscosityModel::New(U.mesh(), phase2Name())),
+    nucModel_(viscosityModel::New(mesh, phase2Name())),
     muModel_(mixtureViscosityModel::New(*this)),
 
-    rhoc_("rho", dimDensity, nucModel_()),
     rhod_("rho", dimDensity, muModel_()),
+    rhoc_("rho", dimDensity, nucModel_()),
 
     alphaMax_(lookupOrDefault("alphaMax", 1.0)),
 
@@ -61,12 +56,12 @@ incompressibleTwoPhaseInteractingMixture
         IOobject
         (
             "rho",
-            U_.time().name(),
-            U_.mesh(),
+            mesh.time().name(),
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        U_.mesh(),
+        mesh,
         dimensionedScalar("rho", dimDensity, 0)
     ),
 
@@ -75,97 +70,37 @@ incompressibleTwoPhaseInteractingMixture
         IOobject
         (
             "nu",
-            U_.time().name(),
-            U_.db()
+            mesh.time().name(),
+            mesh
         ),
-        U_.mesh(),
+        mesh,
         dimensionedScalar(dimViscosity, 0),
         calculatedFvPatchScalarField::typeName
-    )
-{
-    correct();
-}
+    ),
 
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::incompressibleTwoPhaseInteractingMixture::
-~incompressibleTwoPhaseInteractingMixture()
+    Uptr_(nullptr)
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-const Foam::volScalarField&
-Foam::incompressibleTwoPhaseInteractingMixture::alphac() const
-{
-    return alpha2();
-}
-
-
-const Foam::volScalarField&
-Foam::incompressibleTwoPhaseInteractingMixture::alphad() const
-{
-    return alpha1();
-}
-
-
-const Foam::dimensionedScalar&
-Foam::incompressibleTwoPhaseInteractingMixture::rhoc() const
-{
-    return rhoc_;
-};
-
-
-const Foam::dimensionedScalar&
-Foam::incompressibleTwoPhaseInteractingMixture::rhod() const
-{
-    return rhod_;
-}
-
-
-Foam::scalar
-Foam::incompressibleTwoPhaseInteractingMixture::alphaMax() const
-{
-    return alphaMax_;
-}
-
-
-const Foam::volVectorField&
-Foam::incompressibleTwoPhaseInteractingMixture::U() const
-{
-    return U_;
-}
-
-const Foam::volScalarField&
-Foam::incompressibleTwoPhaseInteractingMixture::rho() const
-{
-    return rho_;
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::incompressibleTwoPhaseInteractingMixture::nu() const
-{
-    return nu_;
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::incompressibleTwoPhaseInteractingMixture::nu(const label patchi) const
-{
-    return nu_.boundaryField()[patchi];
-}
-
-
-void Foam::incompressibleTwoPhaseInteractingMixture::correct()
+void Foam::incompressibleDriftFluxMixture::correct()
 {
     rho_ = alpha1()*rhod_ + alpha2()*rhoc_;
-    nu_ = muModel_->mu(rhoc_*nucModel_->nu(), U_)/rho_;
+    nu_ = muModel_->mu(rhoc_*nucModel_->nu(), *Uptr_)/rho_;
 }
 
 
-bool Foam::incompressibleTwoPhaseInteractingMixture::read()
+Foam::incompressibleDriftFluxMixture&
+Foam::incompressibleDriftFluxMixture::initialise(const volVectorField& U)
+{
+    Uptr_ = &U;
+    correct();
+    return *this;
+}
+
+
+bool Foam::incompressibleDriftFluxMixture::read()
 {
     if (twoPhaseMixture::read())
     {

applications/solvers/modules/incompressibleDriftFlux/incompressibleDriftFluxMixture/incompressibleDriftFluxMixture.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,43 +22,40 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::incompressibleTwoPhaseInteractingMixture
+    Foam::incompressibleDriftFluxMixture
 
 Description
-    A two-phase incompressible transportModel for interacting phases
-    requiring the direct evaluation of the mixture viscosity,
-    e.g. activated sludge or slurry.
+    Class to represent a mixture of two constant density phases
 
 SourceFiles
-    incompressibleTwoPhaseInteractingMixture.C
+    incompressibleDriftFluxMixture.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef incompressibleTwoPhaseInteractingMixture_H
-#define incompressibleTwoPhaseInteractingMixture_H
+#ifndef incompressibleDriftFluxMixture_H
+#define incompressibleDriftFluxMixture_H
 
-#include "twoPhaseMixture.H"
+#include "twoPhaseVoFMixture.H"
+#include "incompressibleTwoPhases.H"
 #include "viscosityModel.H"
+#include "mixtureViscosityModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class mixtureViscosityModel;
-
 /*---------------------------------------------------------------------------*\
-          Class incompressibleTwoPhaseInteractingMixture Declaration
+                Class incompressibleDriftFluxMixture Declaration
 \*---------------------------------------------------------------------------*/
 
-class incompressibleTwoPhaseInteractingMixture
+class incompressibleDriftFluxMixture
 :
-    public twoPhaseMixture,
+    public twoPhaseVoFMixture,
+    virtual public incompressibleTwoPhases,
     public viscosity
 {
-    // Private data
-
-        volVectorField& U_;
+    // Private Data
 
         //- Continuous phase viscosity model
         autoPtr<viscosityModel> nucModel_;
@@ -66,12 +63,12 @@ class incompressibleTwoPhaseInteractingMixture
         //- Mixture viscosity model
         autoPtr<mixtureViscosityModel> muModel_;
 
-        //- Continuous phase density
-        dimensionedScalar rhoc_;
-
         //- Dispersed phase density
         dimensionedScalar rhod_;
 
+        //- Continuous phase density
+        dimensionedScalar rhoc_;
+
         //- Optional maximum dispersed phase-fraction (e.g. packing limit)
         scalar alphaMax_;
 
@@ -81,57 +78,100 @@ class incompressibleTwoPhaseInteractingMixture
         //- Mixture kinematic viscosity
         volScalarField nu_;
 
+        //- Pointer to the mixture velocity field
+        //  Set by the initialise(U) function
+        const volVectorField* Uptr_;
+
 
 public:
 
-    TypeName("incompressibleTwoPhaseInteractingMixture");
+    TypeName("incompressibleDriftFluxMixture");
 
 
     // Constructors
 
-        //- Construct from components
-        incompressibleTwoPhaseInteractingMixture
-        (
-            volVectorField& U,
-            const surfaceScalarField& phi
-        );
+        //- Construct from a mesh
+        incompressibleDriftFluxMixture(const fvMesh& mesh);
 
 
     //- Destructor
-    virtual ~incompressibleTwoPhaseInteractingMixture();
+    virtual ~incompressibleDriftFluxMixture()
+    {}
 
 
     // Member Functions
 
-        //- Return const-access to the continuous phase-fraction
-        const volScalarField& alphac() const;
-
         //- Return const-access to the dispersed phase-fraction
-        const volScalarField& alphad() const;
+        const volScalarField& alphad() const
+        {
+            return alpha1();
+        }
 
-        //- Return const-access to continuous-phase density
-        const dimensionedScalar& rhoc() const;
+        //- Return const-access to the continuous phase-fraction
+        const volScalarField& alphac() const
+        {
+            return alpha2();
+        }
 
         //- Return const-access to the dispersed-phase density
-        const dimensionedScalar& rhod() const;
+        const dimensionedScalar& rhod() const
+        {
+            return rhod_;
+        }
+
+        //- Return const-access to continuous-phase density
+        const dimensionedScalar& rhoc() const
+        {
+            return rhoc_;
+        }
+
+        //- Return const-access to the dispersed-phase density
+        const dimensionedScalar& rho1() const
+        {
+            return rhod_;
+        }
+
+        //- Return const-access to continuous-phase density
+        const dimensionedScalar& rho2() const
+        {
+            return rhoc_;
+        }
 
         //- Optional maximum phase-fraction (e.g. packing limit)
         //  Defaults to 1
-        scalar alphaMax() const;
-
-        //- Return const-access to the mixture velocity
-        const volVectorField& U() const;
+        scalar alphaMax() const
+        {
+            return alphaMax_;
+        }
 
         //- Return the mixture density
-        virtual const volScalarField& rho() const;
+        virtual const volScalarField& rho() const
+        {
+            return rho_;
+        }
 
         //- Return the mixture viscosity
-        virtual tmp<volScalarField> nu() const;
+        virtual tmp<volScalarField> nu() const
+        {
+            return nu_;
+        }
 
         //- Return the mixture viscosity for patch
-        virtual tmp<scalarField> nu(const label patchi) const;
+        virtual tmp<scalarField> nu(const label patchi) const
+        {
+            return nu_.boundaryField()[patchi];
+        }
+
+        //- Return the mixture velocity field
+        const volVectorField& U() const
+        {
+            return *Uptr_;
+        }
 
         //- Correct the laminar viscosity
+        incompressibleDriftFluxMixture& initialise(const volVectorField& U);
+
+        //- Correct the density and mixture viscosity
         virtual void correct();
 
         //- Read base phaseProperties dictionary

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "BinghamPlastic.H"
-#include "addToRunTimeSelectionTable.H"
+#include "incompressibleDriftFluxMixture.H"
 #include "fvcGrad.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -49,7 +50,7 @@ namespace mixtureViscosityModels
 
 Foam::mixtureViscosityModels::BinghamPlastic::BinghamPlastic
 (
-    const incompressibleTwoPhaseInteractingMixture& mixture
+    const incompressibleDriftFluxMixture& mixture
 )
 :
     plastic(mixture),

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.H

@@ -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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,7 @@ public:
         //- Construct from mixture
         BinghamPlastic
         (
-            const incompressibleTwoPhaseInteractingMixture& mixture
+            const incompressibleDriftFluxMixture& mixture
         );
 
 

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/Make/files

@@ -5,4 +5,4 @@ BinghamPlastic/BinghamPlastic.C
 slurry/slurry.C
 Quemada/Quemada.C
 
-LIB = $(FOAM_LIBBIN)/libdriftFluxTransportModels
+LIB = $(FOAM_LIBBIN)/libincompressibleDriftFluxViscosityModel

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/Make/options

@@ -1,12 +1,16 @@
 EXE_INC = \
-    -I../incompressibleTwoPhaseInteractingMixture \
+    -I../incompressibleDriftFluxMixture \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhases/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
     -lphysicalProperties \
     -ltwoPhaseMixture \
+    -lincompressibleTwoPhases \
     -lfiniteVolume \
     -lmeshTools

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/Quemada/Quemada.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Quemada.H"
-#include "volFields.H"
+#include "incompressibleDriftFluxMixture.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -49,7 +49,7 @@ namespace mixtureViscosityModels
 
 Foam::mixtureViscosityModels::Quemada::Quemada
 (
-    const incompressibleTwoPhaseInteractingMixture& mixture
+    const incompressibleDriftFluxMixture& mixture
 )
 :
     mixtureViscosityModel(mixture),

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/Quemada/Quemada.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,7 +92,7 @@ public:
     // Constructors
 
         //- Construct from mixture
-        Quemada(const incompressibleTwoPhaseInteractingMixture& mixture);
+        Quemada(const incompressibleDriftFluxMixture& mixture);
 
 
     //- Destructor

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModel.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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "mixtureViscosityModel.H"
+#include "incompressibleDriftFluxMixture.H"
 #include "volFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -39,10 +40,10 @@ namespace Foam
 
 Foam::mixtureViscosityModel::mixtureViscosityModel
 (
-    const incompressibleTwoPhaseInteractingMixture& mixture
+    const incompressibleDriftFluxMixture& mixture
 )
 :
-    viscosityModel(mixture.U().mesh(), mixture.phase1Name()),
+    viscosityModel(mixture.mesh(), mixture.phase1Name()),
     mixture_(mixture)
 {}
 

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModel.H

@@ -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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,6 @@ SourceFiles
 #ifndef mixtureViscosityModel_H
 #define mixtureViscosityModel_H
 
-#include "incompressibleTwoPhaseInteractingMixture.H"
 #include "viscosityModel.H"
 #include "runTimeSelectionTables.H"
 
@@ -51,6 +50,8 @@ SourceFiles
 namespace Foam
 {
 
+class incompressibleDriftFluxMixture;
+
 /*---------------------------------------------------------------------------*\
                     Class mixtureViscosityModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -65,7 +66,7 @@ protected:
     // Protected data
 
         //- Mixture properties
-        const incompressibleTwoPhaseInteractingMixture& mixture_;
+        const incompressibleDriftFluxMixture& mixture_;
 
 
 public:
@@ -82,7 +83,7 @@ public:
             mixtureViscosityModel,
             dictionary,
             (
-                const incompressibleTwoPhaseInteractingMixture& mixture
+                const incompressibleDriftFluxMixture& mixture
             ),
             (mixture)
         );
@@ -93,7 +94,7 @@ public:
         //- Construct from mixture
         mixtureViscosityModel
         (
-            const incompressibleTwoPhaseInteractingMixture& mixture
+            const incompressibleDriftFluxMixture& mixture
         );
 
         //- Disallow default bitwise copy construction
@@ -105,7 +106,7 @@ public:
         //- Return a reference to the selected viscosity model
         static autoPtr<mixtureViscosityModel> New
         (
-            const incompressibleTwoPhaseInteractingMixture& mixture
+            const incompressibleDriftFluxMixture& mixture
         );
 
 

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModelNew.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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,14 +24,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "mixtureViscosityModel.H"
-#include "volFields.H"
-#include "surfaceFields.H"
+#include "incompressibleDriftFluxMixture.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::autoPtr<Foam::mixtureViscosityModel> Foam::mixtureViscosityModel::New
 (
-    const incompressibleTwoPhaseInteractingMixture& mixture
+    const incompressibleDriftFluxMixture& mixture
 )
 {
     const word modelType
@@ -40,7 +39,7 @@ Foam::autoPtr<Foam::mixtureViscosityModel> Foam::mixtureViscosityModel::New
         (
             viscosityModel::findModelDict
             (
-                mixture.U().mesh(),
+                mixture.mesh(),
                 mixture.phase1Name()
             )
         ).lookup("viscosityModel")

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/plastic/plastic.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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "plastic.H"
+#include "incompressibleDriftFluxMixture.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -48,7 +49,7 @@ namespace mixtureViscosityModels
 
 Foam::mixtureViscosityModels::plastic::plastic
 (
-    const incompressibleTwoPhaseInteractingMixture& mixture
+    const incompressibleDriftFluxMixture& mixture
 )
 :
     mixtureViscosityModel(mixture),

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/plastic/plastic.H

@@ -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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,7 +44,7 @@ SourceFiles
 namespace Foam
 {
 
-class incompressibleTwoPhaseInteractingMixture;
+class incompressibleDriftFluxMixture;
 
 namespace mixtureViscosityModels
 {
@@ -83,7 +83,7 @@ public:
     // Constructors
 
         //- Construct from mixture
-        plastic(const incompressibleTwoPhaseInteractingMixture& mixture);
+        plastic(const incompressibleDriftFluxMixture& mixture);
 
 
     //- Destructor

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/slurry/slurry.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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "slurry.H"
+#include "incompressibleDriftFluxMixture.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -48,7 +49,7 @@ namespace mixtureViscosityModels
 
 Foam::mixtureViscosityModels::slurry::slurry
 (
-    const incompressibleTwoPhaseInteractingMixture& mixture
+    const incompressibleDriftFluxMixture& mixture
 )
 :
     mixtureViscosityModel(mixture)

applications/solvers/modules/incompressibleDriftFlux/mixtureViscosityModels/slurry/slurry.H

@@ -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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,7 +53,7 @@ SourceFiles
 namespace Foam
 {
 
-class incompressibleTwoPhaseInteractingMixture;
+class incompressibleDriftFluxMixture;
 
 namespace mixtureViscosityModels
 {
@@ -76,7 +76,7 @@ public:
     // Constructors
 
         //- Construct from mixture
-        slurry(const incompressibleTwoPhaseInteractingMixture& mixture);
+        slurry(const incompressibleDriftFluxMixture& mixture);
 
 
     //- Destructor

applications/solvers/modules/incompressibleDriftFlux/pressureCorrector.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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 "incompressibleDriftFlux.H"
+#include "constrainHbyA.H"
+#include "constrainPressure.H"
+#include "adjustPhi.H"
+#include "findRefCell.H"
+#include "fvcMeshPhi.H"
+#include "fvcFlux.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvcSnGrad.H"
+#include "fvcReconstruct.H"
+#include "fvmLaplacian.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleDriftFlux::pressureCorrector()
+{
+    volVectorField& U = U_;
+    surfaceScalarField& phi(phi_);
+
+    fvVectorMatrix& UEqn = tUEqn.ref();
+    setrAU(UEqn);
+
+    const surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+
+    while (pimple.correct())
+    {
+        volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
+            fvc::flux(HbyA)
+          + 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;
+        }
+    }
+
+    clearrAU();
+    tUEqn.clear();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/Make/files

@@ -3,4 +3,4 @@ simple/simple.C
 general/general.C
 MichaelsBolger/MichaelsBolger.C
 
-LIB = $(FOAM_LIBBIN)/libdriftFluxRelativeVelocityModels
+LIB = $(FOAM_LIBBIN)/libincompressibleDriftFluxRelativeVelocityModels

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/Make/options

@@ -1,13 +1,17 @@
 EXE_INC = \
-    -I../incompressibleTwoPhaseInteractingMixture \
+    -I../incompressibleDriftFluxMixture \
     -I../mixtureViscosityModels/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhases/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
     -lphysicalProperties \
     -ltwoPhaseMixture \
+    -lincompressibleTwoPhases \
     -lfiniteVolume \
     -lmeshTools

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/MichaelsBolger/MichaelsBolger.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,7 @@ namespace relativeVelocityModels
 Foam::relativeVelocityModels::MichaelsBolger::MichaelsBolger
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseInteractingMixture& mixture,
+    const incompressibleDriftFluxMixture& mixture,
     const uniformDimensionedVectorField& g
 )
 :

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/MichaelsBolger/MichaelsBolger.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,7 +96,7 @@ public:
         MichaelsBolger
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseInteractingMixture& mixture,
+            const incompressibleDriftFluxMixture& mixture,
             const uniformDimensionedVectorField& g
         );
 

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/general/general.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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,7 @@ namespace relativeVelocityModels
 Foam::relativeVelocityModels::general::general
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseInteractingMixture& mixture,
+    const incompressibleDriftFluxMixture& mixture,
     const uniformDimensionedVectorField& g
 )
 :

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/general/general.H

@@ -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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,7 +79,7 @@ public:
         general
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseInteractingMixture& mixture,
+            const incompressibleDriftFluxMixture& mixture,
             const uniformDimensionedVectorField& g
         );
 

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C

@@ -72,7 +72,7 @@ Foam::wordList Foam::relativeVelocityModel::UdmPatchFieldTypes() const
 Foam::relativeVelocityModel::relativeVelocityModel
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseInteractingMixture& mixture,
+    const incompressibleDriftFluxMixture& mixture,
     const uniformDimensionedVectorField& g
 )
 :
@@ -83,12 +83,12 @@ Foam::relativeVelocityModel::relativeVelocityModel
         IOobject
         (
             "Udm",
-            mixture.U().time().name(),
-            mixture.U().mesh(),
+            mixture_.time().name(),
+            mixture_.mesh(),
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mixture_.U().mesh(),
+        mixture_.mesh(),
         dimensionedVector(dimVelocity, Zero),
         UdmPatchFieldTypes()
     )
@@ -100,7 +100,7 @@ Foam::relativeVelocityModel::relativeVelocityModel
 Foam::autoPtr<Foam::relativeVelocityModel> Foam::relativeVelocityModel::New
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseInteractingMixture& mixture,
+    const incompressibleDriftFluxMixture& mixture,
     const uniformDimensionedVectorField& g
 )
 {

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.H

@@ -34,7 +34,7 @@ SourceFiles
 #ifndef relativeVelocityModel_H
 #define relativeVelocityModel_H
 
-#include "incompressibleTwoPhaseInteractingMixture.H"
+#include "incompressibleDriftFluxMixture.H"
 #include "uniformDimensionedFields.H"
 #include "runTimeSelectionTables.H"
 
@@ -60,7 +60,7 @@ protected:
     // Protected data
 
         //- Mixture properties
-        const incompressibleTwoPhaseInteractingMixture& mixture_;
+        const incompressibleDriftFluxMixture& mixture_;
 
         //- Acceleration due to gravity
         const uniformDimensionedVectorField& g_;
@@ -82,7 +82,7 @@ public:
         dictionary,
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseInteractingMixture& mixture,
+            const incompressibleDriftFluxMixture& mixture,
             const uniformDimensionedVectorField& g
         ),
         (dict, mixture, g)
@@ -95,7 +95,7 @@ public:
         relativeVelocityModel
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseInteractingMixture& mixture,
+            const incompressibleDriftFluxMixture& mixture,
             const uniformDimensionedVectorField& g
         );
 
@@ -107,7 +107,7 @@ public:
     static autoPtr<relativeVelocityModel> New
     (
         const dictionary& dict,
-        const incompressibleTwoPhaseInteractingMixture& mixture,
+        const incompressibleDriftFluxMixture& mixture,
         const uniformDimensionedVectorField& g
     );
 
@@ -119,7 +119,7 @@ public:
     // Member Functions
 
         //- Mixture properties
-        const incompressibleTwoPhaseInteractingMixture& mixture() const
+        const incompressibleDriftFluxMixture& mixture() const
         {
             return mixture_;
         }

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/simple/simple.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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,7 @@ namespace relativeVelocityModels
 Foam::relativeVelocityModels::simple::simple
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseInteractingMixture& mixture,
+    const incompressibleDriftFluxMixture& mixture,
     const uniformDimensionedVectorField& g
 )
 :

applications/solvers/modules/incompressibleDriftFlux/relativeVelocityModels/simple/simple.H

@@ -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) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -73,7 +73,7 @@ public:
         simple
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseInteractingMixture& mixture,
+            const incompressibleDriftFluxMixture& mixture,
             const uniformDimensionedVectorField& g
         );
 

applications/solvers/modules/twoPhaseVoFSolver/alphaPredictor.C

@@ -33,6 +33,41 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::surfaceScalarField> Foam::solvers::twoPhaseVoFSolver::alphaPhi
+(
+    const surfaceScalarField& phi,
+    const volScalarField& alpha,
+    const dictionary& alphaControls
+)
+{
+    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")
+            }
+        )
+    );
+
+    return fvc::flux
+    (
+        phi,
+        alpha,
+        compressionScheme
+    );
+}
+
+
 void Foam::solvers::twoPhaseVoFSolver::alphaSolve
 (
     const dictionary& alphaControls
@@ -55,25 +90,6 @@ void Foam::solvers::twoPhaseVoFSolver::alphaSolve
         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;
@@ -219,11 +235,11 @@ void Foam::solvers::twoPhaseVoFSolver::alphaSolve
         // Split operator
         tmp<surfaceScalarField> talphaPhi1Un
         (
-            fvc::flux
+            alphaPhi
             (
                 phiCN(),
                 (cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime())(),
-                compressionScheme.rewind()
+                alphaControls
             )
         );
 
@@ -399,8 +415,6 @@ void Foam::solvers::twoPhaseVoFSolver::alphaPredictor()
     {
         alphaSolve(alphaControls);
     }
-
-    mixture.correct();
 }
 
 

applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.C

@@ -155,8 +155,8 @@ void Foam::solvers::twoPhaseVoFSolver::preSolve()
 void Foam::solvers::twoPhaseVoFSolver::prePredictor()
 {
     VoFSolver::prePredictor();
-
     alphaPredictor();
+    mixture.correct();
 }
 
 

applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.H

@@ -119,14 +119,21 @@ private:
         //- Solve for the phase-fractions
         void alphaSolve(const dictionary& alphaControls);
 
-        //- Solve for the phase-fractions
-        void alphaPredictor();
-
 
 protected:
 
     // Protected Member Functions
 
+        virtual tmp<surfaceScalarField> alphaPhi
+        (
+            const surfaceScalarField& phi,
+            const volScalarField& alpha,
+            const dictionary& alphaControls
+        );
+
+        //- Solve for the phase-fractions
+        void alphaPredictor();
+
         //- Adjust the rDeltaT in the vicinity of the interface
         virtual void setInterfaceRDeltaT(volScalarField& rDeltaT);
 

applications/solvers/multiphase/driftFluxFoam/Make/files

@@ -1,4 +0,0 @@
-incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.C
-driftFluxFoam.C
-
-EXE = $(FOAM_APPBIN)/driftFluxFoam

applications/solvers/multiphase/driftFluxFoam/Make/options

@@ -1,27 +0,0 @@
-EXE_INC = \
-    -IincompressibleTwoPhaseInteractingMixture \
-    -ImixtureViscosityModels/lnInclude \
-    -I./relativeVelocityModels/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/physicalProperties/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/VoF \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceCompression/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-    -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude
-
-EXE_LIBS = \
-    -ldriftFluxTransportModels \
-    -ldriftFluxRelativeVelocityModels \
-    -lphysicalProperties \
-    -linterfaceCompression \
-    -ltwoPhaseMixture \
-    -lmomentumTransportModels \
-    -lcompressibleMomentumTransportModels \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsampling \
-    -lfvModels \
-    -lfvConstraints

applications/solvers/multiphase/driftFluxFoam/UEqn.H

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

applications/solvers/multiphase/driftFluxFoam/alphaControls.H

@@ -1,15 +0,0 @@
-const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
-
-label nAlphaCorr(alphaControls.lookup<label>("nAlphaCorr"));
-
-label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
-
-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
-bool alphaApplyPrevCorr
-(
-    alphaControls.lookupOrDefault<Switch>("alphaApplyPrevCorr", false)
-);

applications/solvers/multiphase/driftFluxFoam/alphaEqn.H

@@ -1,122 +0,0 @@
-{
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
-
-    if (MULESCorr)
-    {
-        fvScalarMatrix alpha1Eqn
-        (
-            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
-          + fv::gaussConvectionScheme<scalar>
-            (
-                mesh,
-                phi,
-                upwind<scalar>(mesh, phi)
-            ).fvmDiv(phi, alpha1)
-        );
-
-        solve(alpha1Eqn);
-
-        Info<< "Phase-1 volume fraction = "
-            << alpha1.weightedAverage(mesh.Vsc()).value()
-            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
-            << endl;
-
-        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
-        alphaPhi = talphaPhiUD();
-
-        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
-        {
-            Info<< "Applying the previous iteration correction flux" << endl;
-
-            MULES::correct
-            (
-                geometricOneField(),
-                alpha1,
-                alphaPhi,
-                talphaPhiCorr0.ref(),
-                UniformField<scalar>(mixture.alphaMax()),
-                zeroField()
-            );
-
-            alphaPhi += talphaPhiCorr0();
-        }
-
-        // Cache the upwind-flux
-        talphaPhiCorr0 = talphaPhiUD;
-    }
-
-    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
-    {
-        tmp<surfaceScalarField> talphaPhiUn
-        (
-            fvc::flux
-            (
-                phi,
-                alpha1,
-                alphaScheme
-            )
-          + fvc::flux
-            (
-                phir,
-                alpha1,
-                alpharScheme
-            )
-        );
-
-        if (MULESCorr)
-        {
-            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
-            volScalarField alpha10("alpha10", alpha1);
-
-            MULES::correct
-            (
-                geometricOneField(),
-                alpha1,
-                talphaPhiUn(),
-                talphaPhiCorr.ref(),
-                UniformField<scalar>(mixture.alphaMax()),
-                zeroField()
-            );
-
-            // Under-relax the correction for all but the 1st corrector
-            if (aCorr == 0)
-            {
-                alphaPhi += talphaPhiCorr();
-            }
-            else
-            {
-                alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                alphaPhi += 0.5*talphaPhiCorr();
-            }
-        }
-        else
-        {
-            alphaPhi = talphaPhiUn;
-
-            MULES::explicitSolve
-            (
-                geometricOneField(),
-                alpha1,
-                phi,
-                alphaPhi,
-                UniformField<scalar>(mixture.alphaMax()),
-                zeroField()
-            );
-        }
-    }
-
-    if (alphaApplyPrevCorr && MULESCorr)
-    {
-        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
-    }
-
-    alpha2 = 1.0 - alpha1;
-
-    Info<< "Phase-1 volume fraction = "
-        << alpha1.weightedAverage(mesh.Vsc()).value()
-        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
-        << endl;
-}

applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H

@@ -1,74 +0,0 @@
-#include "alphaControls.H"
-
-{
-    surfaceScalarField alphaPhi
-    (
-        IOobject
-        (
-            "alphaPhi",
-            runTime.name(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(phi.dimensions(), 0)
-    );
-
-    surfaceScalarField phir(fvc::flux(relativeVelocity->Udm()));
-
-    if (nAlphaSubCycles > 1)
-    {
-        dimensionedScalar totalDeltaT = runTime.deltaT();
-        surfaceScalarField alphaPhiSum
-        (
-            IOobject
-            (
-                "alphaPhiSum",
-                runTime.name(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar(phi.dimensions(), 0)
-        );
-
-        for
-        (
-            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-            !(++alphaSubCycle).end();
-        )
-        {
-            #include "alphaEqn.H"
-            alphaPhiSum += (runTime.deltaT()/totalDeltaT)*alphaPhi;
-        }
-
-        alphaPhi = alphaPhiSum;
-    }
-    else
-    {
-        #include "alphaEqn.H"
-    }
-
-    // Apply the diffusion term separately to allow implicit solution
-    // and boundedness of the explicit advection
-    {
-        fvScalarMatrix alpha1Eqn
-        (
-            fvm::ddt(alpha1) - fvc::ddt(alpha1)
-          - fvm::laplacian(turbulence->nut(), alpha1)
-        );
-
-        alpha1Eqn.solve("alpha1Diffusion");
-
-        alphaPhi += alpha1Eqn.flux();
-        alpha2 = 1.0 - alpha1;
-
-        Info<< "Phase-1 volume fraction = "
-            << alpha1.weightedAverage(mesh.Vsc()).value()
-            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
-            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
-            << endl;
-    }
-
-    mixture.correct();
-
-    rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
-}

applications/solvers/multiphase/driftFluxFoam/createFields.H

@@ -1,132 +0,0 @@
-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"
-
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-Info<< "Reading incompressibleTwoPhaseInteractingMixture\n" << endl;
-incompressibleTwoPhaseInteractingMixture mixture(U, phi);
-
-volScalarField& alpha1(mixture.alpha1());
-const volScalarField& rho(mixture.rho());
-
-// Mass flux
-surfaceScalarField rhoPhi
-(
-    IOobject
-    (
-        "rhoPhi",
-        runTime.name(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    fvc::interpolate(rho)*phi
-);
-
-
-// Dispersed phase relative velocity model
-autoPtr<relativeVelocityModel> relativeVelocity
-(
-    relativeVelocityModel::New(mixture, mixture, g)
-);
-
-
-// Construct compressible turbulence model
-autoPtr<compressible::momentumTransportModel> turbulence
-(
-    compressible::momentumTransportModel::New(rho, U, rhoPhi, mixture)
-);
-
-
-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(alpha1.name());
-
-// MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
-
-#include "createMRF.H"
-#include "createFvModels.H"
-#include "createFvConstraints.H"
-
-//- Pointer to the surface momentum field
-//  used for ddtCorr with MRF
-autoPtr<surfaceVectorField> Uf;
-
-if (MRF.size())
-{
-    Info<< "Constructing face momentum Uf" << endl;
-
-    // Ensure the U BCs are up-to-date before constructing Uf
-    U.correctBoundaryConditions();
-
-    Uf = new surfaceVectorField
-    (
-        IOobject
-        (
-            "Uf",
-            runTime.name(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        fvc::interpolate(U)
-    );
-}

applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C

@@ -1,143 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 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
-    driftFluxFoam
-
-Description
-    Solver for 2 incompressible fluids using the mixture approach with the
-    drift-flux approximation for relative motion of the phases.
-
-    Used for simulating the settling of the dispersed phase and other similar
-    separation problems.
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "timeSelector.H"
-#include "CMULES.H"
-#include "subCycle.H"
-#include "incompressibleTwoPhaseInteractingMixture.H"
-#include "relativeVelocityModel.H"
-#include "momentumTransportModel.H"
-#include "compressibleMomentumTransportModels.H"
-#include "pimpleControl.H"
-#include "findRefCell.H"
-#include "pressureReference.H"
-#include "fvModels.H"
-#include "fvConstraints.H"
-
-#include "constrainHbyA.H"
-#include "constrainPressure.H"
-#include "adjustPhi.H"
-
-#include "fvcDdt.H"
-#include "fvcSnGrad.H"
-#include "fvcFlux.H"
-#include "fvcMeshPhi.H"
-#include "fvcReconstruct.H"
-
-#include "fvmDdt.H"
-#include "fvmDiv.H"
-#include "fvmLaplacian.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createTimeControls.H"
-    #include "createFields.H"
-    #include "initContinuityErrs.H"
-
-    volScalarField& alpha2(mixture.alpha2());
-    const dimensionedScalar& rho1 = mixture.rhod();
-    const dimensionedScalar& rho2 = mixture.rhoc();
-
-    turbulence->validate();
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (pimple.run(runTime))
-    {
-        #include "readTimeControls.H"
-        #include "CourantNo.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.userTimeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            mixture.correct();
-
-            fvModels.correct();
-
-            #include "alphaEqnSubCycle.H"
-
-            relativeVelocity->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/driftFluxFoam/gh.H

@@ -1,4 +0,0 @@
-    Info<< "Calculating field g.h\n" << endl;
-    dimensionedScalar ghRef(- mag(g)*hRef);
-    volScalarField gh("gh", (g & mesh.C()) - ghRef);
-    surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef);

applications/solvers/multiphase/driftFluxFoam/pEqn.H

@@ -1,71 +0,0 @@
-{
-    volScalarField rAU("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)
-      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi, Uf)
-    );
-    MRF.makeRelative(phiHbyA);
-    adjustPhi(phiHbyA, U, p_rgh);
-
-    surfaceScalarField phig
-    (
-        (
-          - 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"
-
-    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;
-    }
-
-    // Correct Uf for MRF
-    fvc::correctUf(Uf, U, phi, MRF);
-}

applications/solvers/multiphase/driftFluxFoam/readhRef.H

@@ -1,13 +0,0 @@
-    Info<< "\nReading hRef" << endl;
-    uniformDimensionedScalarField hRef
-    (
-        IOobject
-        (
-            "hRef",
-            runTime.constant(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        dimensionedScalar(dimLength, 0)
-    );

bin/driftFluxFoam

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

tutorials/modules/incompressibleDriftFlux/dahl/constant/phaseProperties

@@ -26,5 +26,6 @@ relativeVelocityModel simple;
     residualAlpha   0;
 }
 
+sigma           0;
 
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/dahl/system/controlDict

@@ -14,7 +14,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     driftFluxFoam;
+application     foamRun;
+
+solver          incompressibleDriftFlux;
 
 startFrom       startTime;
 
@@ -48,6 +50,8 @@ adjustTimeStep  on;
 
 maxCo           5;
 
+maxAlphaCo      1;
+
 maxDeltaT       1;
 
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/dahl/system/fvSchemes

@@ -30,8 +30,9 @@ divSchemes
 
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(tauDm)          Gauss linear;
-    div(phi,alpha) Gauss vanLeer;
-    div(phirb,alpha) Gauss linear;
+
+    div(phi,alpha)      Gauss vanLeer;
+
     div(rhoPhi,k)       Gauss limitedLinear 1;
     div(rhoPhi,epsilon) Gauss limitedLinear 1;
 

tutorials/modules/incompressibleDriftFlux/dahl/system/fvSolution

@@ -36,6 +36,7 @@ solvers
     {
         solver          PCG;
         preconditioner  DIC;
+
         tolerance       1e-6;
         relTol          0;
         minIter         1;
@@ -44,10 +45,10 @@ solvers
     p_rgh
     {
         solver          GAMG;
+        smoother        GaussSeidel;
+
         tolerance       1e-7;
         relTol          0.01;
-        smoother        GaussSeidel;
-        nCellsInCoarsestLevel 20;
     }
 
     p_rghFinal
@@ -56,10 +57,17 @@ solvers
         relTol          0;
     }
 
+    pcorr
+    {
+        $p_rghFinal;
+        tolerance       1e-3;
+    }
+
     "(U|k|epsilon)"
     {
         solver          PBiCGStab;
         preconditioner  DILU;
+
         tolerance       1e-7;
         relTol          0.1;
         minIter         1;

tutorials/modules/incompressibleDriftFlux/mixerVessel2DMRF/constant/phaseProperties

@@ -26,5 +26,6 @@ relativeVelocityModel simple;
     residualAlpha   0;
 }
 
+sigma           0;
 
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/mixerVessel2DMRF/system/controlDict

@@ -14,7 +14,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     driftFluxFoam;
+application     foamRun;
+
+solver          incompressibleDriftFlux;
 
 startFrom       startTime;
 
@@ -48,7 +50,8 @@ adjustTimeStep  yes;
 
 maxCo           1;
 
+maxAlphaCo      1;
+
 maxDeltaT       0.05;
 
-
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/mixerVessel2DMRF/system/fvSchemes

@@ -26,12 +26,13 @@ gradSchemes
 
 divSchemes
 {
-    default         none;
+    default             none;
 
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
     div(tauDm)          Gauss linear;
-    div(phi,alpha) Gauss vanLeer;
-    div(phirb,alpha) Gauss linear;
+
+    div(phi,alpha)      Gauss vanLeer;
+
     div(rhoPhi,k)       Gauss limitedLinear 1;
     div(rhoPhi,epsilon) Gauss limitedLinear 1;
 

tutorials/modules/incompressibleDriftFlux/mixerVessel2DMRF/system/fvSolution

@@ -38,26 +38,19 @@ solvers
     {
         solver          PCG;
         preconditioner  DIC;
+
         tolerance       1e-6;
         relTol          0;
         minIter         1;
     }
 
-    "rho.*"
-    {
-        solver          PCG;
-        preconditioner  DIC;
-        tolerance       1e-8;
-        relTol          0;
-    }
-
     p_rgh
     {
         solver          GAMG;
+        smoother        GaussSeidel;
+
         tolerance       1e-6;
         relTol          0.05;
-        smoother        GaussSeidel;
-        nCellsInCoarsestLevel 20;
     }
 
     p_rghFinal
@@ -66,6 +59,12 @@ solvers
         relTol          0;
     }
 
+    pcorr
+    {
+        $p_rghFinal;
+        tolerance       1e-3;
+    }
+
     "(U|k|epsilon)"
     {
         solver          PBiCGStab;
@@ -91,5 +90,13 @@ PIMPLE
     pRefValue       0;
 }
 
+relaxationFactors
+{
+    equations
+    {
+        ".*"           1;
+    }
+}
+
 
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/tank3D/constant/phaseProperties

@@ -26,4 +26,6 @@ relativeVelocityModel simple;
     residualAlpha   0;
 }
 
+sigma           0;
+
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/tank3D/system/controlDict

@@ -14,7 +14,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     driftFluxFoam;
+application     foamRun;
+
+solver          incompressibleDriftFlux;
 
 startFrom       startTime;