solvers::compressibleVoF: New solver module for compressible two-phase flow with VoF

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

Class
    Foam::solvers::compressibleVoF

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

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

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

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

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

SourceFiles
    compressibleVoF.C

See also
    Foam::solvers::fluidSolver

applications/solvers/modules/fluid/compressibleVoF/Make/files

@@ -0,0 +1,9 @@
+setRDeltaT.C
+moveMesh.C
+alphaPredictor.C
+momentumPredictor.C
+thermophysicalPredictor.C
+pressureCorrector.C
+compressibleVoF.C
+
+LIB = $(FOAM_LIBBIN)/libcompressibleVoF

applications/solvers/modules/fluid/compressibleVoF/Make/options

@@ -1,36 +1,30 @@
 EXE_INC = \
-    -I. \
+    -I$(FOAM_SOLVERS)/modules/fluid/fluidSolver/lnInclude \
     -IcompressibleTwoPhaseMixture \
     -IcompressibleInterPhaseTransportModel/lnInclude \
     -ItwoPhaseChange/lnInclude \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
-    -I$(LIB_SRC)/fvModels/lnInclude
 
-EXE_LIBS = \
+LIB_LIBS = \
+    -lfluidSolver \
     -lcompressibleTwoPhaseMixture \
     -ltwoPhaseSurfaceTension \
     -lfluidThermophysicalModels \
-    -lspecie \
     -ltwoPhaseMixture \
     -ltwoPhaseProperties \
     -linterfaceProperties \
     -lcompressibleTwoPhaseChangeModels \
-    -lmomentumTransportModels \
-    -lcompressibleMomentumTransportModels \
     -lcompressibleInterPhaseTransportModel \
     -lfiniteVolume \
     -lfvModels \
     -lfvConstraints \
-    -lmeshTools
+    -lmeshTools \
+    -lsampling

applications/solvers/modules/fluid/compressibleVoF/VoFSurfaceFilm/VoFPatchTransfer/VoFPatchTransfer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/compressibleVoF/VoFSurfaceFilm/VoFPatchTransfer/VoFPatchTransfer.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/compressibleVoF/alphaPredictor.C

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+#include "interfaceCompression.H"
+#include "CMULES.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleVoF::alphaPredictor()
+{
+    #include "alphaControls.H"
+
+    volScalarField& alpha2(mixture.alpha2());
+
+    const volScalarField& rho1 = mixture.thermo1().rho();
+    const volScalarField& rho2 = mixture.thermo2().rho();
+
+    tmp<surfaceScalarField> talphaPhi1(alphaPhi1);
+
+    if (nAlphaSubCycles > 1)
+    {
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+
+        talphaPhi1 = new surfaceScalarField
+        (
+            IOobject
+            (
+                "alphaPhi1",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(alphaPhi1.dimensions(), 0)
+        );
+
+        surfaceScalarField rhoPhiSum
+        (
+            IOobject
+            (
+                "rhoPhiSum",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(rhoPhi.dimensions(), 0)
+        );
+
+        tmp<volScalarField> trSubDeltaT;
+
+        if (LTS)
+        {
+            trSubDeltaT =
+                fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+        }
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            #include "alphaEqn.H"
+            talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+        }
+
+        alphaPhi1 = talphaPhi1();
+        rhoPhi = rhoPhiSum;
+    }
+    else
+    {
+        #include "alphaEqn.H"
+    }
+
+    contErr =
+        (
+            fvc::ddt(rho)()() + fvc::div(rhoPhi)()()
+          - (fvModels().source(alpha1, rho1)&rho1)()
+          - (fvModels().source(alpha2, rho2)&rho2)()
+        );
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/compressibleVoF/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/compressibleVoF/compressibleTwoPhaseMixture/compressibleTwoPhaseMixture.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/compressibleVoF/compressibleVoF.C

@@ -0,0 +1,368 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+#include "localEulerDdtScheme.H"
+#include "hydrostaticInitialisation.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(compressibleVoF, 0);
+    addToRunTimeSelectionTable(solver, compressibleVoF, fvMesh);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::compressibleVoF::correctCoNum()
+{
+    fluidSolver::correctCoNum(phi);
+
+    alphaCoNum = 0;
+    scalar meanAlphaCoNum = 0;
+
+    if (mesh.nInternalFaces())
+    {
+        const scalarField sumPhi
+        (
+            mixture.nearInterface()().primitiveField()
+            *fvc::surfaceSum(mag(phi))().primitiveField()
+        );
+
+        alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+        meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+    }
+
+    Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+        << " max: " << alphaCoNum << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
+:
+    fluidSolver(mesh),
+
+    U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    ),
+
+    phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(U) & mesh.Sf()
+    ),
+
+    mixture(U, phi),
+
+    alpha1(mixture.alpha1()),
+
+    alphaRestart
+    (
+        typeIOobject<surfaceScalarField>
+        (
+            IOobject::groupName("alphaPhi", alpha1.group()),
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ).headerOk()
+    ),
+
+    buoyancy(mesh),
+
+    p_rgh(buoyancy.p_rgh),
+
+    rho(mixture.rho()),
+
+    dgdt
+    (
+        IOobject
+        (
+            "dgdt",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        alpha1*fvc::div(phi)
+    ),
+
+    pressureReference
+    (
+        mixture.p(),
+        p_rgh,
+        pimple.dict(),
+        false
+    ),
+
+    pMin
+    (
+        "pMin",
+        dimPressure,
+        mixture
+    ),
+
+    rhoPhi
+    (
+        IOobject
+        (
+            "rhoPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    ),
+
+    alphaPhi1
+    (
+        IOobject
+        (
+            IOobject::groupName("alphaPhi", alpha1.group()),
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        phi*fvc::interpolate(alpha1)
+    ),
+
+    K("K", 0.5*magSqr(U)),
+
+    turbulence
+    (
+        rho,
+        U,
+        phi,
+        rhoPhi,
+        alphaPhi1,
+        mixture
+    ),
+
+    phaseChangePtr
+    (
+        compressible::twoPhaseChangeModel::New(mixture)
+    ),
+
+    phaseChange(*phaseChangePtr),
+
+    MRF(mesh)
+{
+    // Read the controls
+    read();
+
+    mesh.schemes().setFluxRequired(p_rgh.name());
+    mesh.schemes().setFluxRequired(alpha1.name());
+
+    if (alphaRestart)
+    {
+        Info << "Restarting alpha" << endl;
+    }
+
+    if (mesh.dynamic())
+    {
+        Info<< "Constructing face momentum Uf" << endl;
+
+        Uf = new surfaceVectorField
+        (
+            IOobject
+            (
+                "Uf",
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            fvc::interpolate(U)
+        );
+    }
+
+    if (transient())
+    {
+        correctCoNum();
+    }
+    else if (LTS)
+    {
+        Info<< "Using LTS" << endl;
+
+        trDeltaT = tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fv::localEulerDdt::rDeltaTName,
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                mesh,
+                dimensionedScalar(dimless/dimTime, 1),
+                extrapolatedCalculatedFvPatchScalarField::typeName
+            )
+        );
+    }
+
+    if (correctPhi)
+    {
+        rAU = new volScalarField
+        (
+            IOobject
+            (
+                "rAU",
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimTime/dimDensity, 1)
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::compressibleVoF::~compressibleVoF()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::solvers::compressibleVoF::maxDeltaT() const
+{
+    const scalar maxAlphaCo
+    (
+        runTime.controlDict().lookup<scalar>("maxAlphaCo")
+    );
+
+    const scalar deltaT = fluidSolver::maxDeltaT();
+
+    if (alphaCoNum > small)
+    {
+        return min
+        (
+            deltaT,
+            maxAlphaCo/(alphaCoNum + small)*runTime.deltaTValue()
+        );
+    }
+    else
+    {
+        return deltaT;
+    }
+}
+
+
+void Foam::solvers::compressibleVoF::prePredictor()
+{
+    fvModels().correct();
+    alphaPredictor();
+    turbulence.correctPhasePhi();
+}
+
+
+void Foam::solvers::compressibleVoF::preSolve()
+{
+    // Read the controls
+    read();
+
+    if (transient())
+    {
+        correctCoNum();
+    }
+    else if (LTS)
+    {
+        setRDeltaT();
+    }
+
+    // Store divU from the previous mesh so that it can be mapped
+    // and used in correctPhi to ensure the corrected phi has the
+    // same divergence
+    if (correctPhi && mesh.topoChanged())
+    {
+        // Construct and register divU for mapping
+        divU = new volScalarField
+        (
+            "divU0",
+            fvc::div(fvc::absolute(phi, U))
+        );
+    }
+
+    fvModels().preUpdateMesh();
+
+    // Update the mesh for topology change, mesh to mesh mapping
+    mesh.update();
+}
+
+
+void Foam::solvers::compressibleVoF::momentumTransportCorrector()
+{
+    if (pimple.transportCorr())
+    {
+        turbulence.correct();
+    }
+}
+
+
+void Foam::solvers::compressibleVoF::thermophysicalTransportCorrector()
+{}
+
+
+void Foam::solvers::compressibleVoF::postSolve()
+{
+    divU.clear();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/compressibleVoF/compressibleVoF.H

@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::solvers::compressibleVoF
+
+Description
+    Solver module for for 2 compressible, non-isothermal immiscible fluids
+    using a VOF (volume of fluid) phase-fraction based interface capturing
+    approach, with optional mesh motion and mesh topology changes including
+    adaptive re-meshing.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Either mixture or two-phase transport modelling may be selected.  In the
+    mixture approach a single laminar, RAS or LES model is selected to model the
+    momentum stress.  In the Euler-Euler two-phase approach separate laminar,
+    RAS or LES selected models are selected for each of the phases.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient and steady simulations.
+
+    Optional fvModels and fvConstraints are provided to enhance the simulation
+    in many ways including adding various sources, Lagrangian
+    particles, surface film etc. and constraining or limiting the solution.
+
+SourceFiles
+    compressibleVoF.C
+
+See also
+    Foam::solvers::fluidSolver
+    Foam::solvers::incompressibleFluid
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleVoF_H
+#define compressibleVoF_H
+
+#include "fluidSolver.H"
+#include "compressibleInterPhaseTransportModel.H"
+#include "twoPhaseChangeModel.H"
+#include "buoyancy.H"
+#include "pressureReference.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class compressibleVoF Declaration
+\*---------------------------------------------------------------------------*/
+
+class compressibleVoF
+:
+    public fluidSolver
+{
+
+protected:
+
+    // Kinematic properties
+
+        //- Velocity field
+        volVectorField U;
+
+        //- Volumetric-flux field
+        surfaceScalarField phi;
+
+
+    // Phase properties
+
+        //- The compressible two-phase mixture
+        compressibleTwoPhaseMixture mixture;
+
+        //- Reference to the primary phase-fraction
+        volScalarField& alpha1;
+
+        //- Switch indicating if this is a restart
+        bool alphaRestart;
+
+        scalar alphaCoNum;
+
+
+    // Thermophysical properties
+
+        //- Buoyancy force
+        solvers::buoyancy buoyancy;
+
+        //- Reference to the buoyant pressure for buoyant cases
+        //  otherwise to the pressure
+        volScalarField& p_rgh;
+
+        //- Reference to the mixture continuity density field
+        const volScalarField& rho;
+
+        //- Compressibility source
+        volScalarField::Internal dgdt;
+
+
+    // Pressure reference
+
+        //- Pressure reference
+        Foam::pressureReference pressureReference;
+
+        //- Minimum pressure
+        dimensionedScalar pMin;
+
+
+    // Kinematic properties
+
+        //- Mass flux field
+        surfaceScalarField rhoPhi;
+
+        // Phase-1 volumetric flux
+        surfaceScalarField alphaPhi1;
+
+        //- Kinetic energy field
+        //  Used in the energy equation
+        volScalarField K;
+
+
+    // Momentum transport
+
+        //- Momentum transport model
+        compressibleInterPhaseTransportModel turbulence;
+
+        autoPtr<compressible::twoPhaseChangeModel> phaseChangePtr;
+
+        compressible::twoPhaseChangeModel& phaseChange;
+
+
+    // Optional models
+
+        //- MRF zone list
+        IOMRFZoneList MRF;
+
+
+    // Cached temporary fields
+
+        tmp<volScalarField> rAU;
+
+        tmp<volScalarField::Internal> contErr;
+
+        //- MULES Correction
+        tmp<surfaceScalarField> talphaPhi1Corr0;
+
+        //- Pointer to the surface momentum field
+        //  used to recreate the flux after mesh-change
+        autoPtr<surfaceVectorField> Uf;
+
+        //- Pointer to the momentum divergence field
+        //  used in correctPhi to ensure the corrected phi has the
+        //  same divergence
+        autoPtr<volScalarField> divU;
+
+        //- Optional LTS reciprocal time-step field
+        tmp<volScalarField> trDeltaT;
+
+        //- Cached momentum matrix
+        //  shared between the momentum predictor and pressure corrector
+        tmp<fvVectorMatrix> tUEqn;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Set rDeltaT for LTS
+        virtual void setRDeltaT();
+
+        //- Correct the cached Courant numbers
+        void correctCoNum();
+
+        //- Solve for the phase-fractions
+        void alphaPredictor();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("compressibleVoF");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        compressibleVoF(fvMesh& mesh);
+
+        //- Disallow default bitwise copy construction
+        compressibleVoF(const compressibleVoF&) = delete;
+
+
+    //- Destructor
+    virtual ~compressibleVoF();
+
+
+    // Member Functions
+
+        //- Return the current maximum time-step for stable solution
+        virtual scalar maxDeltaT() const;
+
+        //- Called at the start of the time-step, before the PIMPLE loop
+        virtual void preSolve();
+
+        //- Called at the start of the PIMPLE loop to move the mesh
+        virtual bool moveMesh();
+
+        //- Called at the start of the PIMPLE loop
+        virtual void prePredictor();
+
+        //- Construct and optionally solve the momentum equation
+        virtual void momentumPredictor();
+
+        //- Construct and solve the energy equation,
+        //  convert to temperature
+        //  and update thermophysical and transport properties
+        virtual void thermophysicalPredictor();
+
+        //- Construct and solve the pressure equation in the PISO loop
+        virtual void pressureCorrector();
+
+        //- Correct the momentum transport modelling
+        //  Newtonian, non-Newtonian or turbulent
+        virtual void momentumTransportCorrector();
+
+        //- Correct the thermophysical transport modelling
+        virtual void thermophysicalTransportCorrector();
+
+        //- Called after the PIMPLE loop at the end of the time-step
+        virtual void postSolve();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const compressibleVoF&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/modules/fluid/compressibleVoF/momentumPredictor.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleVoF::momentumPredictor()
+{
+    tUEqn =
+    (
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - fvm::Sp(contErr(), U)
+      + MRF.DDt(rho, U)
+      + turbulence.divDevTau(U)
+     ==
+        fvModels().source(rho, U)
+    );
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    UEqn.relax();
+
+    fvConstraints().constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - buoyancy.ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        fvConstraints().constrain(U);
+
+        K = 0.5*magSqr(U);
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/compressibleVoF/moveMesh.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+#include "CorrectPhi.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::solvers::compressibleVoF::moveMesh()
+{
+    if (pimple.firstIter() || pimple.moveMeshOuterCorrectors())
+    {
+        if (correctPhi && !divU.valid())
+        {
+            // Construct and register divU for correctPhi
+            divU = new volScalarField
+            (
+                "divU0",
+                fvc::div(fvc::absolute(phi, U))
+            );
+        }
+
+        // Move the mesh
+        mesh.move();
+
+        if (mesh.changing())
+        {
+            buoyancy.moveMesh();
+
+            MRF.update();
+
+            if (correctPhi)
+            {
+                // Calculate absolute flux
+                // from the mapped surface velocity
+                phi = mesh.Sf() & Uf();
+
+                correctUphiBCs(U, phi, true);
+
+                CorrectPhi
+                (
+                    phi,
+                    U,
+                    p_rgh,
+                    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+                    divU(),
+                    pressureReference,
+                    pimple
+                );
+
+                // Make the fluxes relative to the mesh motion
+                fvc::makeRelative(phi, U);
+            }
+
+            mixture.correct();
+
+            meshCourantNo();
+
+            divU.clear();
+
+            return true;
+        }
+
+        divU.clear();
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/compressibleVoF/pressureCorrector.C

@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+#include "constrainHbyA.H"
+#include "constrainPressure.H"
+#include "adjustPhi.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleVoF::pressureCorrector()
+{
+    volScalarField& p = mixture.p();
+
+    const volScalarField& alpha2(mixture.alpha2());
+
+    const volScalarField& rho1 = mixture.rho1();
+    const volScalarField& rho2 = mixture.rho2();
+
+    const volScalarField& psi1 = mixture.thermo1().psi();
+    const volScalarField& psi2 = mixture.thermo2().psi();
+
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    if (rAU.valid())
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+
+    const surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
+
+    while (pimple.correct())
+    {
+        volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
+            fvc::flux(HbyA)
+          + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
+        );
+
+        MRF.makeRelative(phiHbyA);
+
+        surfaceScalarField phig
+        (
+            (
+                mixture.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(phaseChange.Sp_rgh(rho, buoyancy.gh, p_rgh));
+
+        // Make the fluxes relative to the mesh motion
+        fvc::makeRelative(phiHbyA, U);
+
+        tmp<fvScalarMatrix> p_rghEqnComp1;
+        tmp<fvScalarMatrix> p_rghEqnComp2;
+
+        if (pimple.transonic())
+        {
+            const surfaceScalarField rho1f(fvc::interpolate(rho1));
+            const surfaceScalarField rho2f(fvc::interpolate(rho2));
+
+            surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
+            surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
+
+            p_rghEqnComp1 =
+                (
+                    (fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f))/rho1
+                  - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+                  + (alpha1/rho1)
+                   *correction
+                    (
+                        psi1*fvm::ddt(p_rgh)
+                      + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                    )
+                );
+
+            p_rghEqnComp2 =
+                (
+                   (fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f))/rho2
+                 - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+                 + (alpha2/rho2)
+                  *correction
+                   (
+                       psi2*fvm::ddt(p_rgh)
+                     + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+                   )
+               );
+        }
+        else
+        {
+            const surfaceScalarField rho1f(fvc::interpolate(rho1));
+            const surfaceScalarField rho2f(fvc::interpolate(rho2));
+
+            p_rghEqnComp1 =
+                (
+                    (fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f))/rho1
+                  - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+                  + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
+                );
+
+            p_rghEqnComp2 =
+                (
+                   (fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f))/rho2
+                 - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+                 + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
+                );
+        }
+
+        if (mesh.moving())
+        {
+            p_rghEqnComp1.ref() += fvc::div(mesh.phi())*alpha1;
+            p_rghEqnComp2.ref() += fvc::div(mesh.phi())*alpha2;
+        }
+
+        p_rghEqnComp1.ref() *= pos(alpha1);
+        p_rghEqnComp2.ref() *= pos(alpha2);
+
+        p_rghEqnComp1.ref() -=
+            (fvModels().source(alpha1, mixture.thermo1().rho())&rho1)/rho1;
+        p_rghEqnComp2.ref() -=
+            (fvModels().source(alpha2, mixture.thermo2().rho())&rho2)/rho2;
+
+        if (pimple.transonic())
+        {
+            p_rghEqnComp1.ref().relax();
+            p_rghEqnComp2.ref().relax();
+        }
+
+        // Cache p_rgh prior to solve for density update
+        volScalarField p_rgh_0(p_rgh);
+
+        while (pimple.correctNonOrthogonal())
+        {
+            fvScalarMatrix p_rghEqnIncomp
+            (
+                fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
+             == Sp_rgh
+            );
+
+            solve
+            (
+                p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp
+            );
+
+            if (pimple.finalNonOrthogonalIter())
+            {
+                dgdt =
+                (
+                    alpha1*(p_rghEqnComp2 & p_rgh)
+                  - alpha2*(p_rghEqnComp1 & p_rgh)
+                );
+
+                phi = phiHbyA + p_rghEqnIncomp.flux();
+
+                p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*buoyancy.gh, pMin);
+                p_rgh = p - (alpha1*rho1 + alpha2*rho2)*buoyancy.gh;
+                p_rgh.correctBoundaryConditions();
+
+                U = HbyA
+                  + rAU()*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+                U.correctBoundaryConditions();
+                fvConstraints().constrain(U);
+            }
+        }
+
+        // Correct Uf if the mesh is moving
+        fvc::correctUf(Uf, U, fvc::absolute(phi, U), MRF);
+
+        // Update densities from change in p_rgh
+        mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
+        mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
+        mixture.correct();
+
+        // Correct p_rgh for consistency with p and the updated densities
+        p_rgh = p - rho*buoyancy.gh;
+        p_rgh.correctBoundaryConditions();
+    }
+
+    K = 0.5*magSqr(U);
+
+    tUEqn.clear();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/compressibleVoF/setRDeltaT.C

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::solvers::compressibleVoF::setRDeltaT()
+{
+    volScalarField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    const scalar maxCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    const scalar maxAlphaCo
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    const scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    const label nAlphaSpreadIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSpreadIter", 1)
+    );
+
+    const scalar alphaSpreadDiff
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadDiff", 0.2)
+    );
+
+    const scalar alphaSpreadMax
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMax", 0.99)
+    );
+
+    const scalar alphaSpreadMin
+    (
+        pimpleDict.lookupOrDefault<scalar>("alphaSpreadMin", 0.01)
+    );
+
+    const label nAlphaSweepIter
+    (
+        pimpleDict.lookupOrDefault<label>("nAlphaSweepIter", 5)
+    );
+
+    const scalar maxDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great)
+    );
+
+    const scalar minDeltaT
+    (
+        pimpleDict.lookupOrDefault<scalar>("minDeltaT", small)
+    );
+
+    const volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Set the reciprocal time-step from the local Courant number
+    // and maximum and minimum time-steps
+    rDeltaT.ref() = min
+    (
+        1/dimensionedScalar(dimTime, minDeltaT),
+        max
+        (
+            1/dimensionedScalar(dimTime, maxDeltaT),
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V())
+        )
+    );
+
+    if (maxAlphaCo < maxCo)
+    {
+        // Further limit the reciprocal time-step
+        // in the vicinity of the interface
+
+        volScalarField alpha1Bar(fvc::average(alpha1));
+
+        rDeltaT.ref() = max
+        (
+            rDeltaT(),
+            pos0(alpha1Bar() - alphaSpreadMin)
+           *pos0(alphaSpreadMax - alpha1Bar())
+           *fvc::surfaceSum(mag(phi))()()
+           /((2*maxAlphaCo)*mesh.V())
+        );
+    }
+
+    // Update the boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    if (nAlphaSpreadIter > 0)
+    {
+        fvc::spread
+        (
+            rDeltaT,
+            alpha1,
+            nAlphaSpreadIter,
+            alphaSpreadDiff,
+            alphaSpreadMax,
+            alphaSpreadMin
+        );
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
+    }
+
+    Info<< "Smoothed flow time scale min/max = "
+        << gMin(1/rDeltaT.primitiveField())
+        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        pimpleDict.found("rDeltaTDampingCoeff")
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        // Damping coefficient (1-0)
+        const scalar rDeltaTDampingCoeff
+        (
+            pimpleDict.lookup<scalar>("rDeltaTDampingCoeff")
+        );
+
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.primitiveField())
+            << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/compressibleVoF/surfaceTensionModels/liquidProperties/liquidPropertiesSurfaceTension.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/compressibleVoF/surfaceTensionModels/liquidProperties/liquidPropertiesSurfaceTension.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/modules/fluid/compressibleVoF/thermophysicalPredictor.C

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleVoF::thermophysicalPredictor()
+{
+    volScalarField& T = mixture.T();
+    const volScalarField& p = mixture.p();
+    const volScalarField& alpha2(mixture.alpha2());
+
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr(), T)
+      - fvm::laplacian(turbulence.alphaEff(), T)
+      + (
+             mixture.totalInternalEnergy()
+          ?
+            fvc::div(fvc::absolute(phi, U), p)()() // - contErr()/rho*p
+          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
+          - (U()&(fvModels().source(rho, U)&U)()) - contErr()*K
+          :
+            p*fvc::div(fvc::absolute(phi, U))()()
+        )
+       *(
+           alpha1()/mixture.thermo1().Cv()()
+         + alpha2()/mixture.thermo2().Cv()()
+        )
+     ==
+        fvModels().source(rho, T)
+    );
+
+    TEqn.relax();
+
+    fvConstraints().constrain(TEqn);
+
+    TEqn.solve();
+
+    fvConstraints().constrain(T);
+
+    mixture.correctThermo();
+    mixture.correct();
+}
+
+
+// ************************************************************************* //

applications/solvers/modules/fluid/isothermalFluid/Make/files

@@ -1,4 +1,3 @@
-buoyancy/buoyancy.C
 setRDeltaT.C
 moveMesh.C
 prePredictor.C

applications/solvers/multiphase/compressibleInterFoam/Make/files

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

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -1,33 +0,0 @@
-{
-    fvScalarMatrix TEqn
-    (
-        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr, T)
-      - fvm::laplacian(turbulence.alphaEff(), T)
-      + (
-             mixture.totalInternalEnergy()
-          ?
-            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
-          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
-          - (U()&(fvModels.source(rho, U)&U)()) - contErr*K
-          :
-            p*fvc::div(fvc::absolute(phi, U))()()
-        )
-       *(
-           alpha1()/mixture.thermo1().Cv()()
-         + alpha2()/mixture.thermo2().Cv()()
-        )
-     ==
-        fvModels.source(rho, T)
-    );
-
-    TEqn.relax();
-
-    fvConstraints.constrain(TEqn);
-
-    TEqn.solve();
-
-    fvConstraints.constrain(T);
-
-    mixture.correctThermo();
-    mixture.correct();
-}

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

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

applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H

@@ -1,67 +0,0 @@
-tmp<surfaceScalarField> talphaPhi1(alphaPhi1);
-
-if (nAlphaSubCycles > 1)
-{
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-
-    talphaPhi1 = new surfaceScalarField
-    (
-        IOobject
-        (
-            "alphaPhi1",
-            runTime.timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(alphaPhi1.dimensions(), 0)
-    );
-
-    surfaceScalarField rhoPhiSum
-    (
-        IOobject
-        (
-            "rhoPhiSum",
-            runTime.timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(rhoPhi.dimensions(), 0)
-    );
-
-    tmp<volScalarField> trSubDeltaT;
-
-    if (LTS)
-    {
-        trSubDeltaT =
-            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
-    }
-
-    for
-    (
-        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-        !(++alphaSubCycle).end();
-    )
-    {
-        #include "alphaEqn.H"
-        talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
-        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
-    }
-
-    rhoPhi = rhoPhiSum;
-}
-else
-{
-    #include "alphaEqn.H"
-}
-
-const surfaceScalarField& alphaPhi1 = talphaPhi1();
-surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
-
-volScalarField::Internal contErr
-(
-    (
-        fvc::ddt(rho) + fvc::div(rhoPhi)
-      - (fvModels.source(alpha1, rho1)&rho1)
-      - (fvModels.source(alpha2, rho2)&rho2)
-    )()
-);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -1,194 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    compressibleInterFoam
-
-Description
-    Solver for 2 compressible, non-isothermal immiscible fluids using a VOF
-    (volume of fluid) phase-fraction based interface capturing approach,
-    with optional mesh motion and mesh topology changes including adaptive
-    re-meshing.
-
-    The momentum and other fluid properties are of the "mixture" and a single
-    momentum equation is solved.
-
-    Either mixture or two-phase transport modelling may be selected.  In the
-    mixture approach a single laminar, RAS or LES model is selected to model the
-    momentum stress.  In the Euler-Euler two-phase approach separate laminar,
-    RAS or LES selected models are selected for each of the phases.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "interfaceCompression.H"
-#include "CMULES.H"
-#include "EulerDdtScheme.H"
-#include "localEulerDdtScheme.H"
-#include "CrankNicolsonDdtScheme.H"
-#include "subCycle.H"
-#include "compressibleInterPhaseTransportModel.H"
-#include "noPhaseChange.H"
-#include "pimpleControl.H"
-#include "pressureReference.H"
-#include "fvModels.H"
-#include "fvConstraints.H"
-#include "CorrectPhi.H"
-#include "fvcSmooth.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "initContinuityErrs.H"
-    #include "createDyMControls.H"
-    #include "createFields.H"
-    #include "createFieldRefs.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
-    #include "createUfIfPresent.H"
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (pimple.run(runTime))
-    {
-        #include "readDyMControls.H"
-
-        if (LTS)
-        {
-            #include "setRDeltaT.H"
-        }
-        else
-        {
-            #include "CourantNo.H"
-            #include "alphaCourantNo.H"
-            #include "setDeltaT.H"
-        }
-
-        fvModels.preUpdateMesh();
-
-        // Store divU from the previous mesh so that it can be mapped
-        // and used in correctPhi to ensure the corrected phi has the
-        // same divergence
-        tmp<volScalarField> divU;
-
-        if (correctPhi && mesh.topoChanged())
-        {
-            // Construct and register divU for mapping
-            divU = new volScalarField
-            (
-                "divU0",
-                fvc::div(fvc::absolute(phi, U))
-            );
-        }
-
-        // Update the mesh for topology change, mesh to mesh mapping
-        mesh.update();
-
-        runTime++;
-
-        Info<< "Time = " << runTime.userTimeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            if (pimple.firstPimpleIter() || pimple.moveMeshOuterCorrectors())
-            {
-                if (correctPhi && !divU.valid())
-                {
-                    // Construct and register divU for correctPhi
-                    divU = new volScalarField
-                    (
-                        "divU0",
-                        fvc::div(fvc::absolute(phi, U))
-                    );
-                }
-
-                // Move the mesh
-                mesh.move();
-
-                if (mesh.changing())
-                {
-                    gh = (g & mesh.C()) - ghRef;
-                    ghf = (g & mesh.Cf()) - ghRef;
-
-                    MRF.update();
-
-                    if (correctPhi)
-                    {
-                        #include "correctPhi.H"
-                    }
-
-                    mixture.correct();
-
-                    if (checkMeshCourantNo)
-                    {
-                        #include "meshCourantNo.H"
-                    }
-                }
-
-                divU.clear();
-            }
-
-            fvModels.correct();
-
-            #include "alphaControls.H"
-            #include "compressibleAlphaEqnSubCycle.H"
-
-            turbulence.correctPhasePhi();
-
-            #include "UEqn.H"
-            #include "TEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.turbCorr())
-            {
-                turbulence.correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = "
-            << runTime.elapsedCpuTime()
-            << " s\n\n" << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/correctPhi.H

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

applications/solvers/multiphase/compressibleInterFoam/createFieldRefs.H

@@ -1,30 +0,0 @@
-compressible::twoPhaseChangeModel& phaseChange = phaseChangePtr();
-
-volScalarField& alpha2(mixture.alpha2());
-
-const volScalarField& rho1 = mixture.thermo1().rho();
-const volScalarField& rho2 = mixture.thermo2().rho();
-
-volScalarField& p = mixture.p();
-volScalarField& T = mixture.T();
-const volScalarField& psi1 = mixture.thermo1().psi();
-const volScalarField& psi2 = mixture.thermo2().psi();
-
-tmp<volScalarField> rAU;
-
-if (correctPhi)
-{
-    rAU = new volScalarField
-    (
-        IOobject
-        (
-            "rAU",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar(dimTime/dimDensity, 1)
-    );
-}

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -1,110 +0,0 @@
-#include "createRDeltaT.H"
-
-Info<< "Reading field p_rgh\n" << endl;
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "createPhi.H"
-
-Info<< "Constructing compressibleTwoPhaseMixture\n" << endl;
-compressibleTwoPhaseMixture mixture(U, phi);
-
-autoPtr<compressible::twoPhaseChangeModel> phaseChangePtr
-(
-    compressible::twoPhaseChangeModel::New(mixture)
-);
-
-volScalarField& alpha1(mixture.alpha1());
-const volScalarField& rho = mixture.rho();
-
-dimensionedScalar pMin
-(
-    "pMin",
-    dimPressure,
-    mixture
-);
-
-mesh.schemes().setFluxRequired(p_rgh.name());
-mesh.schemes().setFluxRequired(alpha1.name());
-
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-pressureReference pressureReference(mixture.p(), p_rgh, pimple.dict(), false);
-
-
-// Mass flux
-// Initialisation does not matter because rhoPhi is reset after the
-// alpha1 solution before it is used in the U equation.
-surfaceScalarField rhoPhi
-(
-    IOobject
-    (
-        "rhoPhi",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    fvc::interpolate(rho)*phi
-);
-
-volScalarField::Internal dgdt
-(
-    IOobject
-    (
-        "dgdt",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    alpha1*fvc::div(phi)
-);
-
-#include "createAlphaFluxes.H"
-
-// Construct compressible turbulence model
-compressibleInterPhaseTransportModel turbulence
-(
-    rho,
-    U,
-    phi,
-    rhoPhi,
-    alphaPhi1,
-    mixture
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-#include "createMRF.H"
-#include "createFvModels.H"
-#include "createFvConstraints.H"

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -1,164 +0,0 @@
-{
-    if (rAU.valid())
-    {
-        rAU.ref() = 1.0/UEqn.A();
-    }
-    else
-    {
-        rAU = 1.0/UEqn.A();
-    }
-
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
-    volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::flux(HbyA)
-      + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
-    );
-
-    MRF.makeRelative(phiHbyA);
-
-    surfaceScalarField phig
-    (
-        (
-            mixture.surfaceTensionForce()
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
-    );
-
-    phiHbyA += phig;
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
-
-    // Cache the phase change pressure source
-    fvScalarMatrix Sp_rgh(phaseChange.Sp_rgh(rho, gh, p_rgh));
-
-    // Make the fluxes relative to the mesh motion
-    fvc::makeRelative(phiHbyA, U);
-
-    tmp<fvScalarMatrix> p_rghEqnComp1;
-    tmp<fvScalarMatrix> p_rghEqnComp2;
-
-    if (pimple.transonic())
-    {
-        #include "rhofs.H"
-
-        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
-        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
-
-        p_rghEqnComp1 =
-            (
-                (fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f))/rho1
-              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
-              + (alpha1/rho1)
-               *correction
-                (
-                    psi1*fvm::ddt(p_rgh)
-                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
-                )
-            );
-
-        p_rghEqnComp2 =
-            (
-               (fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f))/rho2
-             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
-             + (alpha2/rho2)
-              *correction
-               (
-                   psi2*fvm::ddt(p_rgh)
-                 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
-               )
-           );
-    }
-    else
-    {
-        #include "rhofs.H"
-
-        p_rghEqnComp1 =
-            (
-                (fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f))/rho1
-              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
-              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
-            );
-
-        p_rghEqnComp2 =
-            (
-               (fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f))/rho2
-             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
-             + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
-            );
-    }
-
-    if (mesh.moving())
-    {
-        p_rghEqnComp1.ref() += fvc::div(mesh.phi())*alpha1;
-        p_rghEqnComp2.ref() += fvc::div(mesh.phi())*alpha2;
-    }
-
-    p_rghEqnComp1.ref() *= pos(alpha1);
-    p_rghEqnComp2.ref() *= pos(alpha2);
-
-    p_rghEqnComp1.ref() -=
-        (fvModels.source(alpha1, mixture.thermo1().rho())&rho1)/rho1;
-    p_rghEqnComp2.ref() -=
-        (fvModels.source(alpha2, mixture.thermo2().rho())&rho2)/rho2;
-
-    if (pimple.transonic())
-    {
-        p_rghEqnComp1.ref().relax();
-        p_rghEqnComp2.ref().relax();
-    }
-
-    // Cache p_rgh prior to solve for density update
-    volScalarField p_rgh_0(p_rgh);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix p_rghEqnIncomp
-        (
-            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
-         == Sp_rgh
-        );
-
-        solve
-        (
-            p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp
-        );
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            dgdt =
-            (
-                alpha1*(p_rghEqnComp2 & p_rgh)
-              - alpha2*(p_rghEqnComp1 & p_rgh)
-            );
-
-            phi = phiHbyA + p_rghEqnIncomp.flux();
-
-            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
-            p_rgh.correctBoundaryConditions();
-
-            U = HbyA
-              + rAU()*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
-            U.correctBoundaryConditions();
-            fvConstraints.constrain(U);
-        }
-    }
-
-    // Correct Uf if the mesh is moving
-    fvc::correctUf(Uf, U, fvc::absolute(phi, U), MRF);
-
-    // Update densities from change in p_rgh
-    mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
-    mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
-    mixture.correct();
-
-    // Correct p_rgh for consistency with p and the updated densities
-    p_rgh = p - rho*gh;
-    p_rgh.correctBoundaryConditions();
-
-    K = 0.5*magSqr(U);
-}

bin/compressibleInterFoam

@@ -3,7 +3,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     | Website:  https://openfoam.org
-#   \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 # License
@@ -23,30 +23,24 @@
 #     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 #
 # Script
-#     compressibleInterFilmFoam
+#     compressibleInterFoam
 #
 # Description
-#     Script to inform the user that compressibleInterFilmFoam has been replaced
+#     Script to inform the user that compressibleInterFoam has been superseded
-#     by the more general compressibleInterFoam solver with a special fvModel.
+#     and replaced by the more general compressibleVoF solver module executed by
+#     the foamRun application.
 #
 #------------------------------------------------------------------------------
 
-cat << EOF
+cat <<EOF
 
-The compressibleInterFilmFoam solver has solver has been replaced by the more general
+compressibleInterFoam has been superseded and replaced by the more general
-compressibleInterFoam solver, which now supports surface films using the new
+compressibleVoF solver module executed by the foamRun application:
-VoFSurfaceFilm fvModel.
 
-To run with with surface film create a system/fvModels dictionary
+    foamRun -solver compressibleVoF
-containing the VoFSurfaceFilm specification, e.g.
-
-    VoFSurfaceFilm
-    {
-        type    VoFSurfaceFilm;
-
-        phase   water;
-    }
 
 EOF
 
+foamRun -solver compressibleVoF "$@"
+
 #------------------------------------------------------------------------------

etc/config.sh/bash_completion

@@ -493,31 +493,6 @@ _combinePatchFaces_ ()
 }
 complete -o filenames -o nospace -F _combinePatchFaces_ combinePatchFaces
 
-_compressibleInterFoam_ ()
-{
-    local cur="${COMP_WORDS[COMP_CWORD]}"
-    local prev="${COMP_WORDS[COMP_CWORD-1]}"
-    local line=${COMP_LINE}
-    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
-
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvConstraints -listFvModels -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
-    for o in $used ; do opts="${opts/$o/}" ; done
-    extra=""
-
-    [ "$COMP_CWORD" = 1 ] || \
-    case "$prev" in
-        -case)
-            opts="" ; extra="-d" ;;
-        -fileHandler)
-            opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
-            opts="" ; extra="" ;;
-       *) ;;
-    esac
-    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
-}
-complete -o filenames -o nospace -F _compressibleInterFoam_ compressibleInterFoam
-
 _compressibleMultiphaseInterFoam_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"

src/finiteVolume/Make/files

@@ -458,6 +458,7 @@ $(general)/CorrectPhi/correctUphiBCs.C
 $(general)/pressureReference/pressureReference.C
 $(general)/levelSet/levelSet.C
 $(general)/surfaceToVolVelocity/surfaceToVolVelocity.C
+$(general)/buoyancy/buoyancy.C
 
 solutionControl = $(general)/solutionControl
 $(solutionControl)/solutionControl/solutionControl/solutionControl.C

src/fvConstraints/limitTemperature/limitTemperature.H

@@ -37,7 +37,7 @@ Usage
         selectionMode   all;
 
         // field           T;   // Optional energy/temperature field name
-                                // Set to T for compressibleInterFoam
+                                // Set to T for compressibleVoF
 
         phase           gas; // Optional phase name
 

src/twoPhaseModels/twoPhaseMixture/VoF/setDeltaT.H

@@ -33,15 +33,24 @@ Description
 
 if (adjustTimeStep)
 {
-    scalar deltaT =
+    scalar deltaT = runTime.deltaTValue();
-        min(maxCo/(CoNum + small), maxAlphaCo/(alphaCoNum + small))
+
-       *runTime.deltaTValue();
+    if (CoNum > small)
-    deltaT = min(deltaT, fvModels.maxDeltaT());
+    {
+        deltaT = min(maxCo/CoNum*runTime.deltaTValue(), fvModels.maxDeltaT());
+    }
+
+    if (alphaCoNum > small)
+    {
+        deltaT = min(maxAlphaCo/alphaCoNum*runTime.deltaTValue(), deltaT);
+    }
+
     deltaT = min
     (
         min(deltaT, runTime.deltaTValue() + 0.1*deltaT),
         1.2*runTime.deltaTValue()
     );
+
     runTime.setDeltaT(min(deltaT, maxDeltaT));
 
     Info<< "deltaT = " <<  runTime.deltaTValue() << endl;

tutorials/modules/compressibleVoF/ballValve/system/controlDict

@@ -14,9 +14,11 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     compressibleInterFoam;
+application     foamRun;
 
-startFrom       latestTime;
+solver          compressibleVoF;
+
+startFrom       startTime;
 
 startTime       0;