diff --git a/applications/solvers/modules/Allwmake b/applications/solvers/modules/Allwmake
index 640eaccf22..9b2cb328a7 100755
--- a/applications/solvers/modules/Allwmake
+++ b/applications/solvers/modules/Allwmake
@@ -14,6 +14,8 @@ wmake $targetType VoFSolver
wmake $targetType twoPhaseVoFSolver
incompressibleVoF/Allwmake $targetType $*
compressibleVoF/Allwmake $targetType $*
+wmake $targetType multiphaseVoFSolver
+wmake $targetType incompressibleMultiphaseVoF
multiphaseEuler/Allwmake $targetType $*
wmake $targetType solid
solidDisplacement/Allwmake $targetType $*
diff --git a/applications/solvers/modules/VoFSolver/VoFSolver.H b/applications/solvers/modules/VoFSolver/VoFSolver.H
index fa804ee6d7..22ae747aa0 100644
--- a/applications/solvers/modules/VoFSolver/VoFSolver.H
+++ b/applications/solvers/modules/VoFSolver/VoFSolver.H
@@ -25,18 +25,14 @@ Class
Foam::solvers::VoFSolver
Description
- Solver module base-class for for 2 immiscible fluids using a VOF (volume
- of fluid) phase-fraction based interface capturing approach, with optional
- mesh motion and mesh topology changes including adaptive re-meshing.
+ Base solver module base-class for the solution of immiscible fluids using
+ a VOF (volume of fluid) phase-fraction based interface capturing approach,
+ with optional mesh motion and mesh topology changes including adaptive
+ re-meshing.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
- Either mixture or two-phase transport modelling may be selected. In the
- mixture approach a single laminar, RAS or LES model is selected to model the
- momentum stress. In the Euler-Euler two-phase approach separate laminar,
- RAS or LES selected models are selected for each of the phases.
-
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient and steady simulations.
diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/Make/files b/applications/solvers/modules/incompressibleMultiphaseVoF/Make/files
new file mode 100644
index 0000000000..da6c87f468
--- /dev/null
+++ b/applications/solvers/modules/incompressibleMultiphaseVoF/Make/files
@@ -0,0 +1,7 @@
+incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.C
+incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.C
+alphaPredictor.C
+pressureCorrector.C
+incompressibleMultiphaseVoF.C
+
+LIB = $(FOAM_LIBBIN)/libincompressibleMultiphaseVoF
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/Make/options b/applications/solvers/modules/incompressibleMultiphaseVoF/Make/options
similarity index 56%
rename from applications/solvers/multiphase/multiphaseInterFoam/Make/options
rename to applications/solvers/modules/incompressibleMultiphaseVoF/Make/options
index e120ed62f5..06edbdacba 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/Make/options
+++ b/applications/solvers/modules/incompressibleMultiphaseVoF/Make/options
@@ -1,24 +1,26 @@
EXE_INC = \
- -IincompressibleMultiphaseMixture/lnInclude \
+ -I$(FOAM_SOLVERS)/modules/multiphaseVoFSolver/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)/multiphaseModels/multiphaseProperties/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
+ -I$(LIB_SRC)/MomentumTransportModels/phaseIncompressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
-EXE_LIBS = \
- -lincompressibleMultiphaseMixture \
+LIB_LIBS = \
+ -lmultiphaseVoFSolver \
-lphysicalProperties \
+ -lmultiphaseProperties \
-linterfaceCompression \
- -linterfaceProperties \
-lmomentumTransportModels \
-lincompressibleMomentumTransportModels \
+ -lphaseIncompressibleMomentumTransportModels \
-lfiniteVolume \
+ -lmeshTools \
-lfvModels \
-lfvConstraints \
- -lmeshTools \
-lsampling
diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/alphaPredictor.C b/applications/solvers/modules/incompressibleMultiphaseVoF/alphaPredictor.C
new file mode 100644
index 0000000000..a6a5a18656
--- /dev/null
+++ b/applications/solvers/modules/incompressibleMultiphaseVoF/alphaPredictor.C
@@ -0,0 +1,212 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "incompressibleMultiphaseVoF.H"
+#include "subCycle.H"
+#include "CMULES.H"
+#include "fvcFlux.H"
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+void Foam::solvers::incompressibleMultiphaseVoF::alphaSolve
+(
+ const dictionary& alphaControls
+)
+{
+ const scalar cAlpha(alphaControls.lookup("cAlpha"));
+
+ const word alphaScheme("div(phi,alpha)");
+ const word alpharScheme("div(phirb,alpha)");
+
+ surfaceScalarField phic(mag(phi/mesh.magSf()));
+ phic = min(cAlpha*phic, max(phic));
+
+ UPtrList alphas(phases.size());
+ PtrList alphaPhis(phases.size());
+
+ forAll(phases, phasei)
+ {
+ const incompressibleVoFphase& alpha = phases[phasei];
+
+ alphas.set(phasei, &alpha);
+
+ alphaPhis.set
+ (
+ phasei,
+ new surfaceScalarField
+ (
+ "phi" + alpha.name() + "Corr",
+ fvc::flux
+ (
+ phi,
+ alpha,
+ alphaScheme
+ )
+ )
+ );
+
+ surfaceScalarField& alphaPhi = alphaPhis[phasei];
+
+ forAll(phases, phasej)
+ {
+ incompressibleVoFphase& alpha2 = phases[phasej];
+
+ if (&alpha2 == &alpha) continue;
+
+ surfaceScalarField phir(phic*mixture.nHatf(alpha, alpha2));
+
+ alphaPhi += fvc::flux
+ (
+ -fvc::flux(-phir, alpha2, alpharScheme),
+ alpha,
+ alpharScheme
+ );
+ }
+
+ // Limit alphaPhi for each phase
+ MULES::limit
+ (
+ 1.0/mesh.time().deltaT().value(),
+ geometricOneField(),
+ alpha,
+ phi,
+ alphaPhi,
+ zeroField(),
+ zeroField(),
+ oneField(),
+ zeroField(),
+ false
+ );
+ }
+
+ MULES::limitSum(alphas, alphaPhis, phi);
+
+ rhoPhi = Zero;
+
+ volScalarField sumAlpha
+ (
+ IOobject
+ (
+ "sumAlpha",
+ mesh.time().name(),
+ mesh
+ ),
+ mesh,
+ dimensionedScalar(dimless, 0)
+ );
+
+ forAll(phases, phasei)
+ {
+ incompressibleVoFphase& alpha = phases[phasei];
+ surfaceScalarField& alphaPhi = alphaPhis[phasei];
+
+ MULES::explicitSolve
+ (
+ geometricOneField(),
+ alpha,
+ alphaPhi
+ );
+
+ rhoPhi += alphaPhi*alpha.rho();
+
+ Info<< alpha.name() << " volume fraction, min, max = "
+ << alpha.weightedAverage(mesh.V()).value()
+ << ' ' << min(alpha).value()
+ << ' ' << max(alpha).value()
+ << endl;
+
+ sumAlpha += alpha;
+ }
+
+ Info<< "Phase-sum volume fraction, min, max = "
+ << sumAlpha.weightedAverage(mesh.V()).value()
+ << ' ' << min(sumAlpha).value()
+ << ' ' << max(sumAlpha).value()
+ << endl;
+
+ // Correct the sum of the phase-fractions to avoid 'drift'
+ volScalarField sumCorr(1.0 - sumAlpha);
+ forAll(phases, phasei)
+ {
+ incompressibleVoFphase& alpha = phases[phasei];
+ alpha += alpha*sumCorr;
+ }
+}
+
+
+void Foam::solvers::incompressibleMultiphaseVoF::alphaPredictor()
+{
+ const dictionary& alphaControls = mesh.solution().solverDict("alpha");
+
+ const label nAlphaSubCycles(alphaControls.lookup