solvers::twoPhaseSolver: new base-class for twoPhaseVoFSolver and incompressibleDriftFlux

to separate the interface treatment between VoF and drift-flux and avoid code
duplication.

applications/solvers/modules/Allwmake

@@ -14,6 +14,7 @@ wmake $targetType multicomponentFluid
 wmake $targetType XiFluid
 wmake $targetType shockFluid
 wmake $targetType VoFSolver
+wmake $targetType twoPhaseSolver
 wmake $targetType twoPhaseVoFSolver
 incompressibleVoF/Allwmake $targetType $*
 compressibleVoF/Allwmake $targetType $*

applications/solvers/modules/compressibleVoF/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

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

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/compressibleVoF/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

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

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/compressibleVoF/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

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

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/compressibleVoF/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

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

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/compressibleVoF/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

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

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/compressibleVoF/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

applications/solvers/modules/incompressibleDriftFlux/Make/files

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

applications/solvers/modules/incompressibleDriftFlux/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
@@ -17,7 +17,7 @@ EXE_INC = \
 LIB_LIBS = \
     -lincompressibleDriftFluxViscosityModel \
     -lincompressibleDriftFluxRelativeVelocityModels \
-    -ltwoPhaseVoFSolver \
+    -ltwoPhaseSolver \
     -lphysicalProperties \
     -linterfaceCompression \
     -lincompressibleTwoPhases \

applications/solvers/modules/incompressibleDriftFlux/incompressibleDriftFlux.C

@@ -43,8 +43,23 @@ namespace solvers
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleDriftFlux::correctCoNum()
+{
+    VoFSolver::correctCoNum();
+}
+
+
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
+void Foam::solvers::incompressibleDriftFlux::setInterfaceRDeltaT
+(
+    volScalarField& rDeltaT
+)
+{}
+
+
 void Foam::solvers::incompressibleDriftFlux::correctInterface()
 {}
 
@@ -65,7 +80,7 @@ Foam::solvers::incompressibleDriftFlux::surfaceTensionForce() const
 
 Foam::solvers::incompressibleDriftFlux::incompressibleDriftFlux(fvMesh& mesh)
 :
-    twoPhaseVoFSolver
+    twoPhaseSolver
     (
         mesh,
         autoPtr<twoPhaseVoFMixture>(new incompressibleDriftFluxMixture(mesh))
@@ -73,7 +88,7 @@ Foam::solvers::incompressibleDriftFlux::incompressibleDriftFlux(fvMesh& mesh)
 
     mixture
     (
-        refCast<incompressibleDriftFluxMixture>(twoPhaseVoFSolver::mixture)
+        refCast<incompressibleDriftFluxMixture>(twoPhaseSolver::mixture)
        .initialise(U)
     ),
 
@@ -110,6 +125,11 @@ Foam::solvers::incompressibleDriftFlux::incompressibleDriftFlux(fvMesh& mesh)
     // Read the controls
     readControls();
 
+    if (transient())
+    {
+        correctCoNum();
+    }
+
     if (correctPhi || mesh.topoChanging())
     {
         rAU = new volScalarField
@@ -153,6 +173,12 @@ Foam::solvers::incompressibleDriftFlux::~incompressibleDriftFlux()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::scalar Foam::solvers::incompressibleDriftFlux::maxDeltaT() const
+{
+    return fluidSolver::maxDeltaT();
+}
+
+
 void Foam::solvers::incompressibleDriftFlux::prePredictor()
 {
     VoFSolver::prePredictor();
@@ -196,6 +222,12 @@ void Foam::solvers::incompressibleDriftFlux::prePredictor()
 }
 
 
+void Foam::solvers::incompressibleDriftFlux::pressureCorrector()
+{
+    incompressiblePressureCorrector(p);
+}
+
+
 void Foam::solvers::incompressibleDriftFlux::thermophysicalPredictor()
 {}
 

applications/solvers/modules/incompressibleDriftFlux/incompressibleDriftFlux.H

@@ -46,7 +46,7 @@ SourceFiles
 
 See also
     Foam::solvers::VoFSolver
-    Foam::solvers::twoPhaseVoFSolver
+    Foam::solvers::twoPhaseSolver
     Foam::solvers::compressibleVoF
 
 \*---------------------------------------------------------------------------*/
@@ -54,7 +54,7 @@ See also
 #ifndef incompressibleDriftFlux_H
 #define incompressibleDriftFlux_H
 
-#include "twoPhaseVoFSolver.H"
+#include "twoPhaseSolver.H"
 #include "incompressibleDriftFluxMixture.H"
 #include "relativeVelocityModel.H"
 #include "momentumTransportModel.H"
@@ -73,7 +73,7 @@ namespace solvers
 
 class incompressibleDriftFlux
 :
-    public twoPhaseVoFSolver
+    public twoPhaseSolver
 {
 
 protected:
@@ -105,6 +105,16 @@ protected:
         autoPtr<compressible::momentumTransportModel> momentumTransport;
 
 
+private:
+
+    // Private Member Functions
+
+        //- Correct the cached Courant numbers
+        virtual void correctCoNum();
+
+
+protected:
+
     // Protected Member Functions
 
         //- Return the pressure reference
@@ -140,6 +150,9 @@ protected:
             const dictionary& alphaControls
         );
 
+        //- Adjust the rDeltaT in the vicinity of the interface
+        virtual void setInterfaceRDeltaT(volScalarField& rDeltaT);
+
         //- Calculate the alpha equation sources
         virtual void alphaSuSp
         (
@@ -184,6 +197,9 @@ public:
 
     // Member Functions
 
+        //- Return the current maximum time-step for stable solution
+        virtual scalar maxDeltaT() const;
+
         //- Called at the start of the PIMPLE loop
         virtual void prePredictor();
 

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

@@ -1,6 +1,6 @@
 EXE_INC = \
     -I../incompressibleDriftFluxMixture \
-    -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \

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

@@ -1,7 +1,7 @@
 EXE_INC = \
     -I../incompressibleDriftFluxMixture \
     -I../mixtureViscosityModels/lnInclude \
-    -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \

applications/solvers/modules/incompressibleVoF/Make/files

@@ -1,7 +1,6 @@
 incompressibleTwoPhaseVoFMixture/incompressibleTwoPhaseVoFMixture.C
 incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C
 alphaSuSp.C
-pressureCorrector.C
 incompressibleVoF.C
 
 LIB = $(FOAM_LIBBIN)/libincompressibleVoF

applications/solvers/modules/incompressibleVoF/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

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

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/incompressibleVoF/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

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

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/incompressibleVoF/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \

applications/solvers/modules/incompressibleVoF/incompressibleVoF.C

@@ -147,6 +147,12 @@ void Foam::solvers::incompressibleVoF::prePredictor()
 }
 
 
+void Foam::solvers::incompressibleVoF::pressureCorrector()
+{
+    incompressiblePressureCorrector(p);
+}
+
+
 void Foam::solvers::incompressibleVoF::thermophysicalPredictor()
 {}
 

applications/solvers/modules/incompressibleVoF/pressureCorrector.C

@@ -1,155 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 "incompressibleVoF.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::incompressibleVoF::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
-        (
-            (
-                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/isothermalFilm/fvModels/filmVoFTransfer/Make/options

@@ -9,6 +9,7 @@ EXE_INC = -ggdb3 \
     \
     -I$(FOAM_SOLVERS)/modules/compressibleVoF/lnInclude \
     -I$(FOAM_SOLVERS)/modules/twoPhaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/VoF \

applications/solvers/modules/twoPhaseSolver/Make/files

@@ -0,0 +1,6 @@
+twoPhaseVoFMixture/twoPhaseVoFMixture.C
+alphaPredictor.C
+pressureCorrector.C
+twoPhaseSolver.C
+
+LIB = $(FOAM_LIBBIN)/libtwoPhaseSolver

applications/solvers/modules/twoPhaseSolver/Make/options

@@ -0,0 +1,23 @@
+EXE_INC = \
+    -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)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lVoFSolver \
+    -lphysicalProperties \
+    -linterfaceCompression \
+    -linterfaceProperties \
+    -ltwoPhaseMixture \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvModels \
+    -lfvConstraints \
+    -lsampling

applications/solvers/modules/twoPhaseSolver/alphaPredictor.C

@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "twoPhaseVoFSolver.H"
+#include "twoPhaseSolver.H"
 #include "subCycle.H"
 #include "interfaceCompression.H"
 #include "CMULES.H"
@@ -33,7 +33,7 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::surfaceScalarField> Foam::solvers::twoPhaseVoFSolver::alphaPhi
+Foam::tmp<Foam::surfaceScalarField> Foam::solvers::twoPhaseSolver::alphaPhi
 (
     const surfaceScalarField& phi,
     const volScalarField& alpha,
@@ -68,7 +68,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::solvers::twoPhaseVoFSolver::alphaPhi
 }
 
 
-void Foam::solvers::twoPhaseVoFSolver::alphaSolve
+void Foam::solvers::twoPhaseSolver::alphaSolve
 (
     const dictionary& alphaControls
 )
@@ -364,7 +364,7 @@ void Foam::solvers::twoPhaseVoFSolver::alphaSolve
 }
 
 
-void Foam::solvers::twoPhaseVoFSolver::alphaPredictor()
+void Foam::solvers::twoPhaseSolver::alphaPredictor()
 {
     const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
 

applications/solvers/modules/twoPhaseSolver/pressureCorrector.C

@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "incompressibleDriftFlux.H"
+#include "twoPhaseSolver.H"
 #include "constrainHbyA.H"
 #include "constrainPressure.H"
 #include "adjustPhi.H"
@@ -38,7 +38,10 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::solvers::incompressibleDriftFlux::pressureCorrector()
+void Foam::solvers::twoPhaseSolver::incompressiblePressureCorrector
+(
+    volScalarField& p
+)
 {
     volVectorField& U = U_;
     surfaceScalarField& phi(phi_);

applications/solvers/modules/twoPhaseSolver/twoPhaseSolver.C

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "twoPhaseSolver.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(twoPhaseSolver, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::twoPhaseSolver::twoPhaseSolver
+(
+    fvMesh& mesh,
+    autoPtr<twoPhaseVoFMixture> mixturePtr
+)
+:
+    VoFSolver(mesh, autoPtr<VoFMixture>(mixturePtr.ptr())),
+
+    mixture(refCast<twoPhaseVoFMixture>(VoFSolver::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()
+    ),
+
+    alphaPhi1
+    (
+        IOobject
+        (
+            IOobject::groupName("alphaPhi", alpha1.group()),
+            runTime.name(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        phi*fvc::interpolate(alpha1)
+    )
+{
+    mesh.schemes().setFluxRequired(alpha1.name());
+
+    if (alphaRestart)
+    {
+        Info << "Restarting alpha" << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::twoPhaseSolver::~twoPhaseSolver()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::twoPhaseSolver::preSolve()
+{
+    VoFSolver::preSolve();
+
+    // Do not apply previous time-step mesh compression flux
+    // if the mesh topology changed
+    if (mesh().topoChanged())
+    {
+        talphaPhi1Corr0.clear();
+    }
+}
+
+
+void Foam::solvers::twoPhaseSolver::prePredictor()
+{
+    VoFSolver::prePredictor();
+    alphaPredictor();
+    mixture.correct();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/twoPhaseSolver/twoPhaseSolver.H

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::twoPhaseSolver
+
+Description
+    Solver module base-class for for 2 immiscible fluids, 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
+    twoPhaseSolver.C
+
+See also
+    Foam::solvers::fluidSolver
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef twoPhaseSolver_H
+#define twoPhaseSolver_H
+
+#include "VoFSolver.H"
+#include "twoPhaseVoFMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class twoPhaseSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class twoPhaseSolver
+:
+    public VoFSolver
+{
+
+protected:
+
+    // Phase properties
+
+        //- Reference to the twoPhaseVoFMixture
+        twoPhaseVoFMixture& mixture;
+
+        //- Reference to the phase1-fraction
+        volScalarField& alpha1;
+
+        //- Reference to the phase2-fraction
+        volScalarField& alpha2;
+
+        //- Switch indicating if this is a restart
+        bool alphaRestart;
+
+
+    // Kinematic properties
+
+        // Phase-1 volumetric flux
+        surfaceScalarField alphaPhi1;
+
+
+    // Cached temporary fields
+
+        //- MULES Correction
+        tmp<surfaceScalarField> talphaPhi1Corr0;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Solve for the phase-fractions
+        void alphaSolve(const dictionary& alphaControls);
+
+
+protected:
+
+    // Protected Member Functions
+
+        virtual tmp<surfaceScalarField> alphaPhi
+        (
+            const surfaceScalarField& phi,
+            const volScalarField& alpha,
+            const dictionary& alphaControls
+        );
+
+        //- Solve for the phase-fractions
+        void alphaPredictor();
+
+        //- Calculate the alpha equation sources
+        virtual void alphaSuSp
+        (
+            tmp<volScalarField::Internal>& Su,
+            tmp<volScalarField::Internal>& Sp
+        ) = 0;
+
+        //- Correct the interface properties following mesh-change
+        //  and phase-fraction update
+        virtual void correctInterface() = 0;
+
+        //- Return the interface surface tension force for the momentum equation
+        virtual tmp<surfaceScalarField> surfaceTensionForce() const = 0;
+
+        //- Construct and solve the incompressible pressure equation
+        //  in the PISO loop
+        void incompressiblePressureCorrector(volScalarField& p);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("twoPhaseSolver");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        twoPhaseSolver(fvMesh& mesh, autoPtr<twoPhaseVoFMixture>);
+
+        //- Disallow default bitwise copy construction
+        twoPhaseSolver(const twoPhaseSolver&) = delete;
+
+
+    //- Destructor
+    virtual ~twoPhaseSolver();
+
+
+    // Member Functions
+
+        //- Called at the start of the time-step, before the PIMPLE loop
+        virtual void preSolve();
+
+        //- Called at the start of the PIMPLE loop
+        virtual void prePredictor();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const twoPhaseSolver&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/twoPhaseVoFSolver/Make/files

@@ -1,6 +1,4 @@
-twoPhaseVoFMixture/twoPhaseVoFMixture.C
 setInterfaceRDeltaT.C
-alphaPredictor.C
 twoPhaseVoFSolver.C
 
 LIB = $(FOAM_LIBBIN)/libtwoPhaseVoFSolver

applications/solvers/modules/twoPhaseVoFSolver/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/twoPhaseSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
     -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
@@ -11,8 +12,7 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
-    -lVoFSolver \
+    -ltwoPhaseSolver \
-    -lfluidSolver \
     -lphysicalProperties \
     -linterfaceCompression \
     -linterfaceProperties \

applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.C

@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "twoPhaseVoFSolver.H"
-#include "localEulerDdtScheme.H"
 #include "fvcAverage.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -42,7 +41,7 @@ namespace solvers
 
 void Foam::solvers::twoPhaseVoFSolver::correctCoNum()
 {
-    VoFSolver::correctCoNum();
+    twoPhaseSolver::correctCoNum();
 
     const scalarField sumPhi
     (
@@ -83,47 +82,10 @@ Foam::solvers::twoPhaseVoFSolver::twoPhaseVoFSolver
     autoPtr<twoPhaseVoFMixture> mixturePtr
 )
 :
-    VoFSolver(mesh, autoPtr<VoFMixture>(mixturePtr.ptr())),
+    twoPhaseSolver(mesh, mixturePtr),
 
-    mixture(refCast<twoPhaseVoFMixture>(VoFSolver::mixture_)),
+    interface(mixture, alpha1, alpha2, 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()
-    ),
-
-    interface(mixture, alpha1, alpha2, U),
-
-    alphaPhi1
-    (
-        IOobject
-        (
-            IOobject::groupName("alphaPhi", alpha1.group()),
-            runTime.name(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        phi*fvc::interpolate(alpha1)
-    )
 {
-    mesh.schemes().setFluxRequired(alpha1.name());
-
-    if (alphaRestart)
-    {
-        Info << "Restarting alpha" << endl;
-    }
-
     if (transient())
     {
         correctCoNum();
@@ -137,27 +99,4 @@ Foam::solvers::twoPhaseVoFSolver::~twoPhaseVoFSolver()
 {}
 
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-void Foam::solvers::twoPhaseVoFSolver::preSolve()
-{
-    VoFSolver::preSolve();
-
-    // Do not apply previous time-step mesh compression flux
-    // if the mesh topology changed
-    if (mesh().topoChanged())
-    {
-        talphaPhi1Corr0.clear();
-    }
-}
-
-
-void Foam::solvers::twoPhaseVoFSolver::prePredictor()
-{
-    VoFSolver::prePredictor();
-    alphaPredictor();
-    mixture.correct();
-}
-
-
 // ************************************************************************* //

applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.H

@@ -55,8 +55,7 @@ See also
 #ifndef twoPhaseVoFSolver_H
 #define twoPhaseVoFSolver_H
 
-#include "VoFSolver.H"
+#include "twoPhaseSolver.H"
-#include "twoPhaseVoFMixture.H"
 #include "interfaceProperties.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -72,78 +71,31 @@ namespace solvers
 
 class twoPhaseVoFSolver
 :
-    public VoFSolver
+    public twoPhaseSolver
 {
 
 protected:
 
-    // Phase properties
-
-        //- Reference to the twoPhaseVoFMixture
-        twoPhaseVoFMixture& mixture;
-
-        //- Reference to the phase1-fraction
-        volScalarField& alpha1;
-
-        //- Reference to the phase2-fraction
-        volScalarField& alpha2;
-
-        //- Switch indicating if this is a restart
-        bool alphaRestart;
-
-
     // Interface properties
 
         interfaceProperties interface;
 
 
-    // Kinematic properties
-
-        // Phase-1 volumetric flux
-        surfaceScalarField alphaPhi1;
-
-
-    // Cached temporary fields
-
-        //- MULES Correction
-        tmp<surfaceScalarField> talphaPhi1Corr0;
-
-
 private:
 
     // Private Member Functions
 
         //- Correct the cached Courant numbers
-        void correctCoNum();
+        virtual void correctCoNum();
-
-        //- Solve for the phase-fractions
-        void alphaSolve(const dictionary& alphaControls);
 
 
 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);
 
-        //- Calculate the alpha equation sources
-        virtual void alphaSuSp
-        (
-            tmp<volScalarField::Internal>& Su,
-            tmp<volScalarField::Internal>& Sp
-        ) = 0;
-
         //- Correct the interface properties following mesh-change
         //  and phase-fraction update
         virtual void correctInterface();
@@ -173,12 +125,6 @@ public:
 
     // Member Functions
 
-        //- Called at the start of the time-step, before the PIMPLE loop
-        virtual void preSolve();
-
-        //- Called at the start of the PIMPLE loop
-        virtual void prePredictor();
-
 
     // Member Operators
 

tutorials/modules/incompressibleDriftFlux/dahl/constant/phaseProperties

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

tutorials/modules/incompressibleDriftFlux/dahl/system/controlDict

@@ -50,8 +50,6 @@ adjustTimeStep  on;
 
 maxCo           5;
 
-maxAlphaCo      1;
-
 maxDeltaT       1;
 
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/mixerVessel2DMRF/constant/phaseProperties

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

tutorials/modules/incompressibleDriftFlux/mixerVessel2DMRF/system/controlDict

@@ -50,8 +50,6 @@ adjustTimeStep  yes;
 
 maxCo           1;
 
-maxAlphaCo      1;
-
 maxDeltaT       0.05;
 
 // ************************************************************************* //

tutorials/modules/incompressibleDriftFlux/tank3D/constant/phaseProperties

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