Initial commit

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 OpenFOAM Foundation
+    Copyright (C) 2020-2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "VoFSolidificationMeltingSource.H"
+#include "twoPhaseMixtureThermo.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace fv
+    {
+        defineTypeNameAndDebug(VoFSolidificationMeltingSource, 0);
+
+        addToRunTimeSelectionTable
+        (
+            option,
+            VoFSolidificationMeltingSource,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::fv::VoFSolidificationMeltingSource::update()
+{
+    if (curTimeIndex_ == mesh_.time().timeIndex())
+    {
+        return;
+    }
+
+    if (debug)
+    {
+        Info<< type() << ": " << name_
+            << " - updating solid phase fraction" << endl;
+    }
+
+    alphaSolid_.oldTime();
+
+    const twoPhaseMixtureThermo& thermo
+    (
+        mesh_.lookupObject<twoPhaseMixtureThermo>
+        (
+            twoPhaseMixtureThermo::dictName
+        )
+    );
+
+    const volScalarField& TVoF = thermo.thermo1().T();
+    const volScalarField CpVoF(thermo.thermo1().Cp());
+    const volScalarField& alphaVoF = thermo.alpha1();
+
+    forAll(cells_, i)
+    {
+        const label celli = cells_[i];
+
+        alphaSolid_[celli] = min
+        (
+            relax_*alphaVoF[celli]*alphaSolidT_->value(TVoF[celli])
+          + (1 - relax_)*alphaSolid_[celli],
+            alphaVoF[celli]
+        );
+    }
+
+    alphaSolid_.correctBoundaryConditions();
+
+    curTimeIndex_ = mesh_.time().timeIndex();
+}
+
+
+Foam::word Foam::fv::VoFSolidificationMeltingSource::alphaSolidName() const
+{
+    const twoPhaseMixtureThermo& thermo
+    (
+        mesh_.lookupObject<twoPhaseMixtureThermo>
+        (
+            twoPhaseMixtureThermo::dictName
+        )
+    );
+
+    const volScalarField& alphaVoF = thermo.alpha1();
+
+    return IOobject::groupName(alphaVoF.name(), "solid");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::VoFSolidificationMeltingSource::VoFSolidificationMeltingSource
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    fv::cellSetOption(sourceName, modelType, dict, mesh),
+    alphaSolidT_(Function1<scalar>::New("alphaSolidT", coeffs_, &mesh)),
+    L_("L", dimEnergy/dimMass, coeffs_),
+    relax_(coeffs_.getOrDefault("relax", 0.9)),
+    Cu_(coeffs_.getOrDefault<scalar>("Cu", 100000)),
+    q_(coeffs_.getOrDefault<scalar>("q", 0.001)),
+    alphaSolid_
+    (
+        IOobject
+        (
+            alphaSolidName(),
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(dimless, Zero),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    curTimeIndex_(-1)
+{
+    fieldNames_.resize(2);
+    fieldNames_[0] = "U";
+    fieldNames_[1] = "T";
+
+    fv::option::resetApplied();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    apply(geometricOneField(), eqn);
+}
+
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    apply(rho, eqn);
+}
+
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    fvMatrix<vector>& eqn,
+    const label fieldi
+)
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    update();
+
+    scalarField& Sp = eqn.diag();
+    const scalarField& V = mesh_.V();
+
+    forAll(cells_, i)
+    {
+        const label celli = cells_[i];
+        const scalar Vc = V[celli];
+        const scalar alphaFluid = 1 - alphaSolid_[celli];
+
+        const scalar S = Cu_*sqr(1 - alphaFluid)/(pow3(alphaFluid) + q_);
+
+        Sp[celli] -= Vc*S;
+    }
+}
+
+
+void Foam::fv::VoFSolidificationMeltingSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldi
+)
+{
+    // Momentum source uses a Boussinesq approximation - redirect
+    addSup(eqn, fieldi);
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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::fv::VoFSolidificationMeltingSource
+
+Description
+    Solidification and melting model for VoF simulations.
+
+    The presence of the solid phase in the flow field is incorporated into the
+    model as a momentum porosity contribution; the energy associated with the
+    phase change is added as an enthalpy contribution.  The solid fraction as a
+    function of temperature \c alphaSolidT is specified as a Foam::Function1.
+
+    The model writes the field \c alpha[01].solid which can be visualised to to
+    show the solid distribution.
+
+Usage
+    Example usage:
+    \verbatim
+    VoFSolidificationMeltingSource1
+    {
+        type            VoFSolidificationMeltingSource;
+        active          yes;
+
+        selectionMode   cellZone;
+        cellZone        solidZone;
+
+        alphaSolidT     table
+        (
+            (330 1)
+            (335 0)
+        );
+
+        L               334000;
+    }
+    \endverbatim
+
+    Where:
+    \table
+        Property     | Description             | Required    | Default value
+        alphaSolidT  | Solid fraction as function of temperature | yes |
+        L            | Latent heat of fusion [J/kg] | yes    |
+        relax        | Relaxation coefficient [0-1] | no     | 0.9
+        Cu           | Model coefficient       | no          | 100000
+        q            | Model coefficient       | no          | 0.001
+    \endtable
+
+See also
+    Foam::fv::solidificationMeltingSource
+    Foam::Function1
+
+SourceFiles
+    VoFSolidificationMeltingSource.C
+    VoFSolidificationMeltingSourceIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VoFSolidificationMeltingSource_H
+#define VoFSolidificationMeltingSource_H
+
+#include "fvMesh.H"
+#include "volFields.H"
+#include "cellSetOption.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class VoFSolidificationMeltingSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class VoFSolidificationMeltingSource
+:
+    public fv::cellSetOption
+{
+    // Private data
+
+        //- Solid fraction as a function of temperature
+        autoPtr<Function1<scalar>> alphaSolidT_;
+
+        //- Latent heat of fusion [J/kg]
+        dimensionedScalar L_;
+
+        //- Phase fraction under-relaxation coefficient
+        scalar relax_;
+
+        //- Mushy region momentum sink coefficient [1/s]; default = 10^5
+        scalar Cu_;
+
+        //- Coefficient used in porosity calc - default = 0.001
+        scalar q_;
+
+        //- Solid phase fraction
+        volScalarField alphaSolid_;
+
+        //- Current time index (used for updating)
+        label curTimeIndex_;
+
+
+    // Private Member Functions
+
+        //- Return the name of the solid phase fraction
+        word alphaSolidName() const;
+
+        //- Update the model
+        void update();
+
+        //- Helper function to apply to the energy equation
+        template<class RhoFieldType>
+        void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn);
+
+        //- No copy construct
+        VoFSolidificationMeltingSource
+        (
+            const VoFSolidificationMeltingSource&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const VoFSolidificationMeltingSource&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("VoFSolidificationMeltingSource");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        VoFSolidificationMeltingSource
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Member Functions
+
+        // Add explicit and implicit contributions
+
+            //- Add explicit contribution to enthalpy equation
+            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldi);
+
+            //- Add implicit contribution to momentum equation
+            virtual void addSup(fvMatrix<vector>& eqn, const label fieldi);
+
+
+        // Add explicit and implicit contributions to compressible equation
+
+            //- Add explicit contribution to compressible enthalpy equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn,
+                const label fieldi
+            );
+
+            //- Add implicit contribution to compressible momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldi
+            );
+
+
+        // IO
+
+            //- Read source dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "VoFSolidificationMeltingSourceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "VoFSolidificationMeltingSource.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::fv::VoFSolidificationMeltingSource::read(const dictionary& dict)
+{
+    if (fv::cellSetOption::read(dict))
+    {
+        alphaSolidT_ = Function1<scalar>::New("alphaSolidT", coeffs_, &mesh_);
+        coeffs_.readEntry("L", L_);
+        coeffs_.readIfPresent("relax", relax_);
+        coeffs_.readIfPresent("Cu", Cu_);
+        coeffs_.readIfPresent("q", q_);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceTemplates.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "fvcDdt.H"
+#include "twoPhaseMixtureThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::fv::VoFSolidificationMeltingSource::apply
+(
+    const RhoFieldType& rho,
+    fvMatrix<scalar>& eqn
+)
+{
+    if (debug)
+    {
+        Info<< type() << ": applying source to " << eqn.psi().name() << endl;
+    }
+
+    update();
+
+    const twoPhaseMixtureThermo& thermo
+    (
+        mesh_.lookupObject<twoPhaseMixtureThermo>
+        (
+            twoPhaseMixtureThermo::dictName
+        )
+    );
+
+    const volScalarField CpVoF(thermo.thermo1().Cp());
+
+    if (eqn.psi().dimensions() == dimTemperature)
+    {
+        eqn += L_/CpVoF*(fvc::ddt(rho, alphaSolid_));
+    }
+    else
+    {
+        eqn += L_*(fvc::ddt(rho, alphaSolid_));
+    }
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/compressibleInterFilmFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 two compressible, non-isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+#include "compressibleInterPhaseTransportModel.H"
+#include "pimpleControl.H"
+#include "SLGThermo.H"
+#include "surfaceFilmModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Solver for two compressible, non-isothermal immiscible fluids"
+        " using VOF phase-fraction based interface capturing."
+    );
+
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createFields.H"
+    #include "createSurfaceFilmModel.H"
+
+    volScalarField& p = mixture.p();
+    volScalarField& T = mixture.T();
+    const volScalarField& psi1 = mixture.thermo1().psi();
+    const volScalarField& psi2 = mixture.thermo2().psi();
+
+    regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
+
+    if (!LTS)
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNo.H"
+            #include "alphaCourantNo.H"
+            #include "setDeltaT.H"
+        }
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        surfaceFilm.evolve();
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "compressibleAlphaEqnSubCycle.H"
+
+            turbulence.correctPhasePhi();
+
+            volScalarField::Internal Srho(surfaceFilm.Srho());
+            contErr -= posPart(Srho);
+
+            #include "UEqn.H"
+            #include "TEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence.correct();
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/createSurfaceFilmModel.H

+Info<< "\nConstructing surface film model" << endl;
+
+SLGThermo slgThermo(mesh, mixture.thermo1());
+autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
+(
+    regionModels::surfaceFilmModel::New(mesh, g)
+);

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H

+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi))
+    );
+    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);
+
+    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 =
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/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_rghEqnComp1.ref().relax();
+
+        p_rghEqnComp2 =
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/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)
+               )
+           );
+        p_rghEqnComp2.ref().relax();
+    }
+    else
+    {
+        #include "rhofs.H"
+
+        p_rghEqnComp1 =
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
+            ) - surfaceFilm.Srho()/rho1;
+
+        p_rghEqnComp2 =
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/rho2
+             - fvc::ddt(alpha2) - fvc::div(alphaPhi2)
+             + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh))
+            );
+    }
+
+    // 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)
+        );
+
+        solve
+        (
+            p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+
+            dgdt =
+            (
+                alpha1*(p_rghEqnComp2 & p_rgh)
+              - alpha2*(p_rghEqnComp1 & p_rgh)
+            );
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    // 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));
+
+    rho = alpha1*rho1 + alpha2*rho2;
+
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
+    K = 0.5*magSqr(U);
+}

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) OpenCFD OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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
+
+Group
+    grpMultiphaseSolvers
+
+Description
+    Solver for two compressible, non-isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+
+    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 "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+#include "compressibleInterPhaseTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Solver for two compressible, non-isothermal immiscible fluids"
+        " using VOF phase-fraction based interface capturing."
+    );
+
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createFields.H"
+
+    volScalarField& p = mixture.p();
+    volScalarField& T = mixture.T();
+    const volScalarField& psi1 = mixture.thermo1().psi();
+    const volScalarField& psi2 = mixture.thermo2().psi();
+
+    if (!LTS)
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNo.H"
+            #include "alphaCourantNo.H"
+            #include "setDeltaT.H"
+        }
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "compressibleAlphaEqnSubCycle.H"
+
+            turbulence.correctPhasePhi();
+
+            #include "UEqn.H"
+            volScalarField divUp("divUp", fvc::div(fvc::absolute(phi, U), p));
+            #include "TEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence.correct();
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/Make/files

+compressibleInterIsoFoam.C
+
+EXE = $(FOAM_APPBIN)/compressibleInterIsoFoam

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/Make/options

+EXE_INC = \
+    -I.. \
+    -I../../VoF \
+    -I../twoPhaseMixtureThermo \
+    -I../VoFphaseCompressibleTurbulenceModels \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/geometricVoF/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lsurfMesh \
+    -lmeshTools \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltwoPhaseMixtureThermo \
+    -ltwoPhaseSurfaceTension \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -ltwoPhaseMixture \
+    -ltwoPhaseProperties \
+    -linterfaceProperties \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lVoFphaseCompressibleTurbulenceModels \
+    -lgeometricVoF

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaControls.H

+const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
+label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaEqn.H

+// Update alpha1
+#include "alphaSuSp.H"
+advector.advect(Sp,(Su + divU*min(alpha1(), scalar(1)))());
+
+// Update rhoPhi
+rhoPhi = advector.getRhoPhi(rho1, rho2);
+alphaPhi10 = advector.alphaPhi();
+
+alpha2 = 1.0 - alpha1;
+
+Info<< "Phase-1 volume fraction = "
+    << alpha1.weightedAverage(mesh.Vsc()).value()
+    << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+    << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
+    << endl;

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaSuSp.H

+volScalarField::Internal Sp
+(
+    IOobject
+    (
+        "Sp",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(dgdt.dimensions(), Zero)
+);
+
+volScalarField::Internal Su
+(
+    IOobject
+    (
+        "Su",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(dgdt.dimensions(), Zero)
+);
+
+forAll(dgdt, celli)
+{
+    if (dgdt[celli] > 0.0)
+    {
+        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+    }
+    else if (dgdt[celli] < 0.0)
+    {
+        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
+    }
+}
+
+volScalarField::Internal divU
+(
+    mesh.moving()
+  ? fvc::div(phi + mesh.phi())
+  : fvc::div(phi)
+);

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/compressibleAlphaEqnSubCycle.H

+if (pimple.nCorrPIMPLE() > 1)
+{
+    if (!pimple.firstIter())
+    {
+        // Resetting alpha1 to value before advection in first PIMPLE
+        // iteration.
+        alpha1 = alpha1.oldTime();
+    }
+}
+tmp<surfaceScalarField> talphaPhi1(alphaPhi10);
+
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+
+    talphaPhi1 = new surfaceScalarField
+    (
+        IOobject
+        (
+            "alphaPhi1",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar(alphaPhi10.dimensions(), Zero)
+    );
+
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar(rhoPhi.dimensions(), Zero)
+    );
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi10;
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + alpha2*rho2;
+
+const surfaceScalarField& alphaPhi1 = talphaPhi1();
+surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
+
+volScalarField::Internal contErr
+(
+    (
+        fvc::ddt(rho) + fvc::div(rhoPhi)
+      - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+      - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+    )()
+);

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/compressibleInterIsoFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020 Johan Roenby
+    Copyright (C) 2020 DLR
+-------------------------------------------------------------------------------
+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
+    compressibleInterFlow
+
+Description
+    Solver derived from interFoam for two compressible, immiscible
+    fluids using the isoAdvector phase-fraction based interface capturing
+    approach, with optional mesh motion and mesh topology changes including
+    adaptive re-meshing.
+
+    Reference:
+    \verbatim
+        Roenby, J., Bredmose, H. and Jasak, H. (2016).
+        A computational method for sharp interface advection
+        Royal Society Open Science, 3
+        doi 10.1098/rsos.160405
+    \endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "CMULES.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+#include "compressibleInterPhaseTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+#include "fvcSmooth.H"
+#include "dynamicRefineFvMesh.H"
+#include "isoAdvection.H"
+#include "twoPhaseMixtureThermo.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Solver for two compressible, non-isothermal immiscible fluids"
+        " using VOF phase-fraction based interface capturing.\n"
+        "With optional mesh motion and mesh topology changes including"
+        " adaptive re-meshing."
+    );
+
+    #include "postProcess.H"
+
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "createUf.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    volScalarField& p = mixture.p();
+    volScalarField& T = mixture.T();
+    const volScalarField& psi1 = mixture.thermo1().psi();
+    const volScalarField& psi2 = mixture.thermo2().psi();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+
+        // Store divU and divUp from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
+
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+
+        ++runTime;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+                if (isA<dynamicRefineFvMesh>(mesh))
+                {
+                    advector.surf().reconstruct();
+                }
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+
+                    if (isA<dynamicRefineFvMesh>(mesh))
+                    {
+                        advector.surf().mapAlphaField();
+                        alpha2 = 1.0 - alpha1;
+                        alpha2.correctBoundaryConditions();
+                        rho == alpha1*rho1 + alpha2*rho2;
+                        rho.correctBoundaryConditions();
+                        rho.oldTime() = rho;
+                        alpha2.oldTime() = alpha2;
+                    }
+
+                    MRF.update();
+
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
+
+                }
+
+                if ((mesh.changing() && correctPhi))
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
+
+                    #include "correctPhi.H"
+
+                    // Make the fluxes relative to the mesh motion
+                    fvc::makeRelative(phi, U);
+
+                    mixture.correct();
+                }
+
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
+            }
+
+            #include "alphaControls.H"
+            #include "compressibleAlphaEqnSubCycle.H"
+
+            turbulence.correctPhasePhi();
+
+            #include "UEqn.H"
+            volScalarField divUp("divUp", fvc::div(fvc::absolute(phi, U), p));
+            #include "TEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence.correct();
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/correctPhi.H

+CorrectPhi
+(
+    U,
+    phi,
+    p,
+    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    divU,
+    pimple
+);
+
+//***HGW phi.oldTime() = phi;
+
+#include "continuityErrs.H"

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/createFields.H

+#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 twoPhaseMixtureThermo\n" << endl;
+twoPhaseMixtureThermo mixture(U, phi);
+
+
+volScalarField& alpha1(mixture.alpha1());
+volScalarField& alpha2(mixture.alpha2());
+
+Info<< "Reading thermophysical properties\n" << endl;
+
+const volScalarField& rho1 = mixture.thermo1().rho();
+const volScalarField& rho2 = mixture.thermo2().rho();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    alpha1*rho1 + alpha2*rho2
+);
+
+
+dimensionedScalar pMin
+(
+    "pMin",
+    dimPressure,
+    mixture
+);
+
+mesh.setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(alpha1.name());
+
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+// 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 dgdt(alpha1*fvc::div(phi));
+
+#include "createAlphaFluxes.H"
+
+Foam::isoAdvection advector(alpha1,phi,U);
+// Construct compressible turbulence model
+compressibleInterPhaseTransportModel turbulence
+(
+    rho,
+    U,
+    phi,
+    rhoPhi,
+    alphaPhi10,
+    mixture
+);
+
+#include "createK.H"
+
+#include "createMRF.H"
+#include "createFvOptions.H"

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/pEqn.H

+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, 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);
+
+    // 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)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/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_rghEqnComp1.ref().relax();
+
+        p_rghEqnComp2 =
+            (
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/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)
+               )
+           );
+        p_rghEqnComp2.ref().relax();
+    }
+    else
+    {
+        #include "rhofs.H"
+
+        p_rghEqnComp1 =
+            pos(alpha1)
+           *(
+                (
+                    fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1*rho1f)
+                  - (fvOptions(alpha1, mixture.thermo1().rho())&rho1)
+                )/rho1
+              - fvc::ddt(alpha1) - fvc::div(alphaPhi1)
+              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh))
+            );
+
+        p_rghEqnComp2 =
+            pos(alpha2)
+           *(
+               (
+                   fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2*rho2f)
+                 - (fvOptions(alpha2, mixture.thermo2().rho())&rho2)
+               )/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);
+
+    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)
+        );
+
+         solve
+         (
+             p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
+             mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+         );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+
+            dgdt =
+            (
+                alpha1*(p_rghEqnComp2 & p_rgh)
+              - alpha2*(p_rghEqnComp1 & p_rgh)
+            );
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    // Correct Uf if the mesh is moving
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += n*(fvc::absolute(phi, U)/mesh.magSf() - (n & Uf));
+    }
+
+
+    // 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));
+
+    rho = alpha1*rho1 + alpha2*rho2;
+
+    // Correct p_rgh for consistency with p and the updated densities
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+
+    K = 0.5*magSqr(U);
+}

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/createFields.H

+#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 twoPhaseMixtureThermo\n" << endl;
+twoPhaseMixtureThermo mixture(U, phi);
+
+volScalarField& alpha1(mixture.alpha1());
+volScalarField& alpha2(mixture.alpha2());
+
+Info<< "Reading thermophysical properties\n" << endl;
+
+const volScalarField& rho1 = mixture.thermo1().rho();
+const volScalarField& rho2 = mixture.thermo2().rho();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    alpha1*rho1 + alpha2*rho2
+);
+
+
+dimensionedScalar pMin
+(
+    "pMin",
+    dimPressure,
+    mixture
+);
+
+mesh.setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(alpha1.name());
+
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+// 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 dgdt(alpha1*fvc::div(phi));
+
+#include "createAlphaFluxes.H"
+
+// Construct compressible turbulence model
+compressibleInterPhaseTransportModel turbulence
+(
+    rho,
+    U,
+    phi,
+    rhoPhi,
+    alphaPhi10,
+    mixture
+);
+
+#include "createK.H"
+
+#include "createMRF.H"
+#include "createFvOptions.H"

OpenFOAM-v2112/applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/Make/files

+overCompressibleInterDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/overCompressibleInterDyMFoam