solvers::incompressibleMultiphaseVoF: New solver module for multiphase VoF simulations

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

incompressibleMultiphaseVoF is derived from the multiphaseVoFSolver which adds
multiphase capability to the VoFSolver base-class used as the basis of all
two-phase and multiphase VoF solvers.

Class
    Foam::solvers::incompressibleMultiphaseVoF

Description
    Solver module for the solution of multiple incompressible, isothermal
    immiscible fluids using a VOF (volume of fluid) phase-fraction based
    interface capturing approach, with optional mesh motion and mesh topology
    changes including adaptive re-meshing.

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

    A mixture approach for momentum transport is provided 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.

SourceFiles
    incompressibleMultiphaseVoF.C

See also
    Foam::solvers::VoFSolver
    Foam::solvers::multiphaseVoFSolver

applications/solvers/modules/multiphaseVoFSolver/Make/files

@@ -0,0 +1,7 @@
+multiphaseVoFMixture/VoFphase/VoFphase.C
+multiphaseVoFMixture/multiphaseVoFMixture.C
+setInterfaceRDeltaT.C
+#alphaPredictor.C
+multiphaseVoFSolver.C
+
+LIB = $(FOAM_LIBBIN)/libmultiphaseVoFSolver

applications/solvers/modules/multiphaseVoFSolver/Make/options

@@ -0,0 +1,20 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
+    -I$(LIB_SRC)/physicalProperties/lnInclude \
+    -I$(LIB_SRC)/multiphaseModels/multiphaseProperties/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lVoFSolver \
+    -lfluidSolver \
+    -lphysicalProperties \
+    -lmultiphaseProperties \
+    -linterfaceCompression \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvModels \
+    -lfvConstraints \
+    -lsampling

applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFMixture/VoFphase/VoFphase.C

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "VoFphase.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::VoFphase::VoFphase(const word& phaseName, const fvMesh& mesh)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phaseName),
+            mesh.time().name(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    ),
+    name_(phaseName)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::VoFphase> Foam::VoFphase::clone() const
+{
+    NotImplemented;
+    return autoPtr<VoFphase>(nullptr);
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFMixture/VoFphase/VoFphase.H

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::VoFphase
+
+Description
+    Single incompressible VoFphase derived from the phase-fraction.
+    Used as part of the multiPhaseMixture for interface-capturing multi-VoFphase
+    simulations.
+
+SourceFiles
+    VoFphase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VoFphase_H
+#define VoFphase_H
+
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class VoFphase Declaration
+\*---------------------------------------------------------------------------*/
+
+class VoFphase
+:
+    public volScalarField
+{
+    // Private Data
+
+        word name_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        VoFphase
+        (
+            const word& name,
+            const fvMesh& mesh
+        );
+
+        //- Return clone
+        virtual autoPtr<VoFphase> clone() const;
+
+        //- Return a pointer to a new VoFphase
+        //  created on freestore
+        //  from Istream
+        class iNew
+        {
+        protected:
+
+            const fvMesh& mesh_;
+
+        public:
+
+            iNew(const fvMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            virtual ~iNew()
+            {}
+
+            virtual autoPtr<VoFphase> operator()(Istream& is) const = 0;
+        };
+
+
+    // Member Functions
+
+        const word& name() const
+        {
+            return name_;
+        }
+
+        const word& keyword() const
+        {
+            return name();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFMixture/multiphaseVoFMixture.C

@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "multiphaseVoFMixture.H"
+#include "correctContactAngle.H"
+#include "surfaceInterpolate.H"
+#include "fvcGrad.H"
+#include "fvcSnGrad.H"
+#include "fvcDiv.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(multiphaseVoFMixture, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::multiphaseVoFMixture::calcAlphas()
+{
+    scalar level = 0.0;
+    alphas_ == 0.0;
+
+    forAll(phases_, phasei)
+    {
+        alphas_ += level*phases_[phasei];
+        level += 1.0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseVoFMixture::multiphaseVoFMixture
+(
+    const fvMesh& mesh,
+    const VoFphase::iNew& inew
+)
+:
+    VoFMixture(mesh),
+    IOdictionary
+    (
+        IOobject
+        (
+            "phaseProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+
+    phases_(lookup("phases"), inew),
+
+    mesh_(mesh),
+
+    alphas_
+    (
+        IOobject
+        (
+            "alphas",
+            mesh_.time().name(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar(dimless, 0)
+    ),
+
+    sigmas_(lookup("sigmas")),
+    dimSigma_(1, 0, -2, 0, 0),
+
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh_.V()), 1.0/3.0)
+    )
+{
+    calcAlphas();
+    alphas_.write();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseVoFMixture::surfaceTensionForce(const volVectorField& U) const
+{
+    tmp<surfaceScalarField> tstf
+    (
+        surfaceScalarField::New
+        (
+            "surfaceTensionForce",
+            mesh_,
+            dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), 0)
+        )
+    );
+
+    surfaceScalarField& stf = tstf.ref();
+
+    forAll(phases_, phasei)
+    {
+        const VoFphase& alpha1 = phases_[phasei];
+
+        for (label phasej = phasei+1; phasej<phases_.size(); phasej++)
+        {
+            const VoFphase& alpha2 = phases_[phasej];
+
+            sigmaTable::const_iterator sigma =
+                sigmas_.find(interfacePair(alpha1, alpha2));
+
+            if (sigma == sigmas_.end())
+            {
+                FatalErrorInFunction
+                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << " in list of sigma values"
+                    << exit(FatalError);
+            }
+
+            stf += dimensionedScalar(dimSigma_, sigma())
+               *fvc::interpolate(K(alpha1, alpha2, U))*
+                (
+                    fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
+                  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
+                );
+        }
+    }
+
+    return tstf;
+}
+
+
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseVoFMixture::nHatfv
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    /*
+    // Cell gradient of alpha
+    volVectorField gradAlpha =
+        alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
+    */
+
+    surfaceVectorField gradAlphaf
+    (
+        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
+      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
+    );
+
+    // Face unit interface normal
+    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseVoFMixture::nHatf
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    // Face unit interface normal flux
+    return nHatfv(alpha1, alpha2) & mesh_.Sf();
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseVoFMixture::K
+(
+    const VoFphase& alpha1,
+    const VoFphase& alpha2,
+    const volVectorField& U
+) const
+{
+    tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
+
+    correctContactAngle
+    (
+        alpha1,
+        alpha2,
+        U.boundaryField(),
+        deltaN_,
+        tnHatfv.ref().boundaryFieldRef()
+    );
+
+    // Simple expression for curvature
+    return -fvc::div(tnHatfv & mesh_.Sf());
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseVoFMixture::nearInterface() const
+{
+    tmp<volScalarField> tnearInt
+    (
+        volScalarField::New
+        (
+            "nearInterface",
+            mesh_,
+            dimensionedScalar(dimless, 0)
+        )
+    );
+
+    forAll(phases_, phasei)
+    {
+        tnearInt.ref() = max
+        (
+            tnearInt(),
+            pos0(phases_[phasei] - 0.01)*pos0(0.99 - phases_[phasei])
+        );
+    }
+
+    return tnearInt;
+}
+
+
+bool Foam::multiphaseVoFMixture::read()
+{
+    if (regIOobject::read())
+    {
+        lookup("sigmas") >> sigmas_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFMixture/multiphaseVoFMixture.H

@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::multiphaseVoFMixture
+
+Description
+    Multiphase VoF mixture with support for interface properties
+
+    Surface tension and contact-angle is handled for the interface
+    between each phase-pair.
+
+SourceFiles
+    multiphaseVoFMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseVoFMixture_H
+#define multiphaseVoFMixture_H
+
+#include "VoFMixture.H"
+#include "IOdictionary.H"
+#include "PtrListDictionary.H"
+#include "VoFphase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class multiphaseVoFMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseVoFMixture
+:
+    public VoFMixture,
+    public IOdictionary
+{
+public:
+
+    class interfacePair
+    :
+        public Pair<word>
+    {
+    public:
+
+        class hash
+        :
+            public Hash<interfacePair>
+        {
+        public:
+
+            hash()
+            {}
+
+            label operator()(const interfacePair& key) const
+            {
+                return word::hash()(key.first()) + word::hash()(key.second());
+            }
+        };
+
+
+        // Constructors
+
+            interfacePair()
+            {}
+
+            interfacePair(const word& alpha1Name, const word& alpha2Name)
+            :
+                Pair<word>(alpha1Name, alpha2Name)
+            {}
+
+            interfacePair(const VoFphase& alpha1, const VoFphase& alpha2)
+            :
+                Pair<word>(alpha1.name(), alpha2.name())
+            {}
+
+
+        // Friend Operators
+
+            friend bool operator==
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return
+                (
+                    ((a.first() == b.first()) && (a.second() == b.second()))
+                 || ((a.first() == b.second()) && (a.second() == b.first()))
+                );
+            }
+
+            friend bool operator!=
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return (!(a == b));
+            }
+    };
+
+
+protected:
+
+    // Protected Data
+
+        //- Dictionary of phases
+        PtrListDictionary<VoFphase> phases_;
+
+        const fvMesh& mesh_;
+
+        volScalarField alphas_;
+
+        typedef HashTable<scalar, interfacePair, interfacePair::hash>
+            sigmaTable;
+
+        sigmaTable sigmas_;
+        dimensionSet dimSigma_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+
+    // Private Member Functions
+
+        void calcAlphas();
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<volScalarField> K
+        (
+            const VoFphase& alpha1,
+            const VoFphase& alpha2,
+            const volVectorField& U
+        ) const;
+
+
+public:
+
+    TypeName("multiphaseVoFMixture");
+
+
+    // Constructors
+
+        //- Construct from fvMesh
+        multiphaseVoFMixture
+        (
+            const fvMesh& mesh,
+            const VoFphase::iNew&
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseVoFMixture()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        PtrListDictionary<VoFphase>& phases()
+        {
+            return phases_;
+        }
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> surfaceTensionForce
+        (
+            const volVectorField& U
+        ) const;
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Read base phaseProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFSolver.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "multiphaseVoFSolver.H"
+#include "localEulerDdtScheme.H"
+#include "CorrectPhi.H"
+#include "geometricZeroField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(multiphaseVoFSolver, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::multiphaseVoFSolver::correctCoNum()
+{
+    VoFSolver::correctCoNum();
+
+    const scalarField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    const scalar meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+
+    Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+        << " max: " << alphaCoNum << endl;
+}
+
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::solvers::multiphaseVoFSolver::correctInterface()
+{}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::solvers::multiphaseVoFSolver::surfaceTensionForce() const
+{
+    return mixture.surfaceTensionForce(U);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::multiphaseVoFSolver::multiphaseVoFSolver
+(
+    fvMesh& mesh,
+    autoPtr<multiphaseVoFMixture> mixturePtr
+)
+:
+    VoFSolver(mesh, autoPtr<VoFMixture>(mixturePtr.ptr())),
+
+    mixture(refCast<multiphaseVoFMixture>(VoFSolver::mixture)),
+
+    phases(mixture.phases())
+{
+    if (transient())
+    {
+        correctCoNum();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::multiphaseVoFSolver::~multiphaseVoFSolver()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::multiphaseVoFSolver::preSolve()
+{
+    VoFSolver::preSolve();
+}
+
+
+void Foam::solvers::multiphaseVoFSolver::prePredictor()
+{
+    VoFSolver::prePredictor();
+
+    alphaPredictor();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFSolver.H

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::multiphaseVoFSolver
+
+Description
+    Base solver module for the solution of multiple immiscible fluids using a
+    VOF (volume of fluid) phase-fraction based interface capturing approach,
+    with optional mesh motion and mesh topology changes including adaptive
+    re-meshing.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient and steady simulations.
+
+SourceFiles
+    multiphaseVoFSolver.C
+
+See also
+    Foam::solvers::VoFSolver
+    Foam::solvers::incompressibleMultiphaseVoFSolver
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseVoFSolver_H
+#define multiphaseVoFSolver_H
+
+#include "VoFSolver.H"
+#include "multiphaseVoFMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class multiphaseVoFSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseVoFSolver
+:
+    public VoFSolver
+{
+
+protected:
+
+    // Phase properties
+
+        //- Reference to the multiphaseVoFMixture
+        multiphaseVoFMixture& mixture;
+
+        //- Reference to the phases
+        PtrListDictionary<VoFphase>& phases;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Correct the cached Courant numbers
+        void correctCoNum();
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Adjust the rDeltaT in the vicinity of the interface
+        virtual void setInterfaceRDeltaT(volScalarField& rDeltaT);
+
+        //- Solve for the phase-fractions
+        virtual void alphaPredictor() = 0;
+
+        //- Correct the interface properties following mesh-change
+        //  and phase-fraction update
+        virtual void correctInterface();
+
+        //- Return the interface surface tension force for the momentum equation
+        virtual tmp<surfaceScalarField> surfaceTensionForce() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseVoFSolver");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        multiphaseVoFSolver(fvMesh& mesh, autoPtr<multiphaseVoFMixture>);
+
+        //- Disallow default bitwise copy construction
+        multiphaseVoFSolver(const multiphaseVoFSolver&) = delete;
+
+
+    //- Destructor
+    virtual ~multiphaseVoFSolver();
+
+
+    // 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 multiphaseVoFSolver&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/multiphaseVoFSolver/setInterfaceRDeltaT.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "multiphaseVoFSolver.H"
+#include "fvcSmooth.H"
+#include "fvcSurfaceIntegrate.H"
+#include "fvcAverage.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::multiphaseVoFSolver::setInterfaceRDeltaT
+(
+    volScalarField& rDeltaT
+)
+{
+    const dictionary& pimpleDict = pimple.dict();
+
+    const scalar maxCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    const scalar maxAlphaCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    const scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    const scalar alphaSpreadMax
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
+    );
+
+    const scalar alphaSpreadMin
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
+    );
+
+    if (maxAlphaCo < maxCo)
+    {
+        // Further limit the reciprocal time-step
+        // in the vicinity of the interface
+
+        const volScalarField::Internal alphaCoRdeltaT
+        (
+            fvc::surfaceSum(mag(phi))()()/((2*maxAlphaCo)*mesh.V())
+        );
+
+        forAll(phases, phasei)
+        {
+            const volScalarField alphaBar(fvc::average(phases[phasei]));
+
+            rDeltaT.ref() = max
+            (
+                rDeltaT(),
+                pos0(alphaBar() - alphaSpreadMin)
+               *pos0(alphaSpreadMax - alphaBar())
+               *alphaCoRdeltaT
+            );
+        }
+    }
+
+    // Update the boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "Flow and interface time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    Info<< "Smoothed flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+}
+
+
+// ************************************************************************* //