diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/Make/files b/applications/solvers/multiphase/compressibleLesInterFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..0ec38b8d7a0c420f902dd2d25e2ab0d109c17d95
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/Make/files
@@ -0,0 +1,3 @@
+compressibleLesInterFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/compressibleLesInterFoam
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/Make/options b/applications/solvers/multiphase/compressibleLesInterFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..7f7317e77a1e286b8cb5bc39c371044c9b592126
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/Make/options
@@ -0,0 +1,15 @@
+INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
+
+EXE_INC = \
+ -I$(LIB_SRC)/transportModels \
+ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+ -I$(LIB_SRC)/LESmodels \
+ -I$(LIB_SRC)/LESmodels/LESdeltas/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+ -linterfaceProperties \
+ -lincompressibleTransportModels \
+ -lincompressibleLESmodels \
+ -lfiniteVolume
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..089591736773862c21899dfc34679e0166748b83
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H
@@ -0,0 +1,29 @@
+ surfaceScalarField muf =
+ twoPhaseProperties.muf()
+ + fvc::interpolate(rho*turbulence->nuSgs());
+
+ fvVectorMatrix UEqn
+ (
+ fvm::ddt(rho, U)
+ + fvm::div(rhoPhi, U)
+ - fvm::laplacian(muf, U)
+ - (fvc::grad(U) & fvc::grad(muf))
+ //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
+ );
+
+ if (momentumPredictor)
+ {
+ solve
+ (
+ UEqn
+ ==
+ fvc::reconstruct
+ (
+ (
+ fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+ - ghf*fvc::snGrad(rho)
+ - fvc::snGrad(pd)
+ ) * mesh.magSf()
+ )
+ );
+ }
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqns.H b/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqns.H
new file mode 100644
index 0000000000000000000000000000000000000000..12eea31127e74f8a38b9ce7f0eab662ff3b1d9f4
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqns.H
@@ -0,0 +1,76 @@
+{
+ word alphaScheme("div(phi,alpha)");
+ word alpharScheme("div(phirb,alpha)");
+
+ surfaceScalarField phir = phic*interface.nHatf();
+
+ for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+ {
+ volScalarField::DimensionedInternalField Sp
+ (
+ IOobject
+ (
+ "Sp",
+ runTime.timeName(),
+ mesh
+ ),
+ mesh,
+ dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
+ );
+
+ volScalarField::DimensionedInternalField Su
+ (
+ IOobject
+ (
+ "Su",
+ runTime.timeName(),
+ mesh
+ ),
+ // Divergence term is handled explicitly to be
+ // consistent with the explicit transport solution
+ divU*min(alpha1, 1.0)
+ );
+
+ forAll(dgdt, celli)
+ {
+ if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
+ {
+ Sp[celli] -= dgdt[celli]*alpha1[celli];
+ Su[celli] += dgdt[celli]*alpha1[celli];
+ }
+ else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
+ {
+ Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
+ }
+ }
+
+
+ surfaceScalarField phiAlpha1 =
+ fvc::flux
+ (
+ phi,
+ alpha1,
+ alphaScheme
+ )
+ + fvc::flux
+ (
+ -fvc::flux(-phir, alpha2, alpharScheme),
+ alpha1,
+ alpharScheme
+ );
+
+ MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
+
+ surfaceScalarField rho1f = fvc::interpolate(rho1);
+ surfaceScalarField rho2f = fvc::interpolate(rho2);
+ rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
+
+ alpha2 = scalar(1) - alpha1;
+ }
+
+ Info<< "Liquid phase volume fraction = "
+ << alpha1.weightedAverage(mesh.V()).value()
+ << " Min(alpha1) = " << min(alpha1).value()
+ << " Min(alpha2) = " << min(alpha2).value()
+ << endl;
+}
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqnsSubCycle.H
new file mode 100644
index 0000000000000000000000000000000000000000..c52dce96900965eefc1571d18fda1aa6ed714c84
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqnsSubCycle.H
@@ -0,0 +1,43 @@
+{
+ label nAlphaCorr
+ (
+ readLabel(piso.lookup("nAlphaCorr"))
+ );
+
+ label nAlphaSubCycles
+ (
+ readLabel(piso.lookup("nAlphaSubCycles"))
+ );
+
+ surfaceScalarField phic = mag(phi/mesh.magSf());
+ phic = min(interface.cGamma()*phic, max(phic));
+
+ volScalarField divU = fvc::div(phi);
+
+ if (nAlphaSubCycles > 1)
+ {
+ dimensionedScalar totalDeltaT = runTime.deltaT();
+ surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+
+ for
+ (
+ subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+ !(++alphaSubCycle).end();
+ )
+ {
+# include "alphaEqns.H"
+ rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+ }
+
+ rhoPhi = rhoPhiSum;
+ }
+ else
+ {
+# include "alphaEqns.H"
+ }
+
+ if (oCorr == 0)
+ {
+ interface.correct();
+ }
+}
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/compressibleLesInterFoam.C b/applications/solvers/multiphase/compressibleLesInterFoam/compressibleLesInterFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..6a7874f1968fab4e166c394b5fcfad0aa6adbafc
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/compressibleLesInterFoam.C
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
+ \\/ 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 2 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, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+ compressibleLesInterFoam
+
+Description
+ Solver for 2 compressible, 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 is modelled using a run-time
+ selectable incompressible LES model.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "interfaceProperties.H"
+#include "twoPhaseMixture.H"
+#include "incompressible/LESmodel/LESmodel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ #include "setRootCase.H"
+ #include "createTime.H"
+ #include "createMesh.H"
+ #include "readEnvironmentalProperties.H"
+ #include "readControls.H"
+ #include "initContinuityErrs.H"
+ #include "createFields.H"
+ #include "CourantNo.H"
+ #include "setInitialDeltaT.H"
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ while (runTime.run())
+ {
+ #include "readControls.H"
+ #include "CourantNo.H"
+ #include "setDeltaT.H"
+
+ runTime++;
+
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ turbulence->correct();
+
+ // --- Outer-corrector loop
+ for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+ {
+ #include "alphaEqnsSubCycle.H"
+
+ solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+
+ #include "UEqn.H"
+
+ // --- PISO loop
+ for (int corr=0; corr<nCorr; corr++)
+ {
+ #include "pEqn.H"
+ }
+ }
+
+ rho = alpha1*rho1 + alpha2*rho2;
+
+ runTime.write();
+
+ Info<< "ExecutionTime = "
+ << runTime.elapsedCpuTime()
+ << " s\n\n" << endl;
+ }
+
+ Info<< "End\n" << endl;
+
+ return(0);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/createFields.H b/applications/solvers/multiphase/compressibleLesInterFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..75260de6c26ff95acd7205e439134857cc2f51d0
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/createFields.H
@@ -0,0 +1,152 @@
+ Info<< "Reading field pd\n" << endl;
+ volScalarField pd
+ (
+ IOobject
+ (
+ "pd",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+ Info<< "Reading field alpha1\n" << endl;
+ volScalarField alpha1
+ (
+ IOobject
+ (
+ "alpha1",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+ Info<< "Calculating field alpha1\n" << endl;
+ volScalarField alpha2("alpha2", scalar(1) - alpha1);
+
+ Info<< "Reading field U\n" << endl;
+ volVectorField U
+ (
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+ #include "createPhi.H"
+
+
+ Info<< "Calculating field g.h\n" << endl;
+ volScalarField gh("gh", g & mesh.C());
+ surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+
+ Info<< "Reading transportProperties\n" << endl;
+ twoPhaseMixture twoPhaseProperties(U, phi);
+
+ dimensionedScalar rho10
+ (
+ twoPhaseProperties.subDict
+ (
+ twoPhaseProperties.phase1Name()
+ ).lookup("rho0")
+ );
+
+ dimensionedScalar rho20
+ (
+ twoPhaseProperties.subDict
+ (
+ twoPhaseProperties.phase2Name()
+ ).lookup("rho0")
+ );
+
+ dimensionedScalar psi1
+ (
+ twoPhaseProperties.subDict
+ (
+ twoPhaseProperties.phase1Name()
+ ).lookup("psi")
+ );
+
+ dimensionedScalar psi2
+ (
+ twoPhaseProperties.subDict
+ (
+ twoPhaseProperties.phase2Name()
+ ).lookup("psi")
+ );
+
+ dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
+
+ volScalarField p
+ (
+ IOobject
+ (
+ "p",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ max
+ (
+ (pd + gh*(alpha1*rho10 + alpha2*rho20))
+ /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
+ pMin
+ )
+ );
+
+ volScalarField rho1 = rho10 + psi1*p;
+ volScalarField rho2 = rho20 + psi2*p;
+
+ volScalarField rho
+ (
+ IOobject
+ (
+ "rho",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ alpha1*rho1 + alpha2*rho2
+ );
+
+
+ // 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
+ (
+ "rho*phi",
+ runTime.timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ fvc::interpolate(rho)*phi
+ );
+
+ volScalarField dgdt =
+ pos(alpha2)*fvc::div(phi)/max(alpha2, 0.0001);
+
+ // Construct interface from alpha1 distribution
+ interfaceProperties interface(alpha1, U, twoPhaseProperties);
+
+ // Construct LES model
+ autoPtr<LESmodel> turbulence
+ (
+ LESmodel::New(U, phi, twoPhaseProperties)
+ );
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..b794d4df2d57dfeb3e78e2eeef2016055a080d7a
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H
@@ -0,0 +1,71 @@
+{
+ volScalarField rUA = 1.0/UEqn.A();
+ surfaceScalarField rUAf = fvc::interpolate(rUA);
+
+ tmp<fvScalarMatrix> pdEqnComp;
+
+ if (transonic)
+ {
+ pdEqnComp =
+ (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+ }
+ else
+ {
+ pdEqnComp =
+ (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+ }
+
+
+ U = rUA*UEqn.H();
+
+ surfaceScalarField phiU
+ (
+ "phiU",
+ (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
+ );
+
+ phi = phiU +
+ (
+ fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+ - ghf*fvc::snGrad(rho)
+ )*rUAf*mesh.magSf();
+
+ for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+ {
+ fvScalarMatrix pdEqnIncomp
+ (
+ fvc::div(phi)
+ - fvm::laplacian(rUAf, pd)
+ );
+
+ solve
+ (
+ (max(alpha1, 0.0)*(psi1/rho1) + max(alpha2, 0.0)*(psi2/rho2))
+ *pdEqnComp()
+ + pdEqnIncomp
+ );
+
+ if (nonOrth == nNonOrthCorr)
+ {
+ dgdt =
+ (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
+ *(pdEqnComp & pd);
+ phi += pdEqnIncomp.flux();
+ }
+ }
+
+ U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+ U.correctBoundaryConditions();
+
+ p = max
+ (
+ (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
+ pMin
+ );
+
+ rho1 = rho10 + psi1*p;
+ rho2 = rho20 + psi2*p;
+
+ Info<< "max(U) " << max(mag(U)).value() << endl;
+ Info<< "min(pd) " << min(pd).value() << endl;
+}
diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/readControls.H b/applications/solvers/multiphase/compressibleLesInterFoam/readControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..7e23354f47faf454cde356329b7bae37d39b2cd9
--- /dev/null
+++ b/applications/solvers/multiphase/compressibleLesInterFoam/readControls.H
@@ -0,0 +1,20 @@
+ #include "readPISOControls.H"
+ #include "readTimeControls.H"
+
+ label nAlphaCorr
+ (
+ readLabel(piso.lookup("nAlphaCorr"))
+ );
+
+ label nAlphaSubCycles
+ (
+ readLabel(piso.lookup("nAlphaSubCycles"))
+ );
+
+ if (nAlphaSubCycles > 1 && nOuterCorr != 1)
+ {
+ FatalErrorIn(args.executable())
+ << "Sub-cycling alpha is only allowed for PISO, "
+ "i.e. when the number of outer-correctors = 1"
+ << exit(FatalError);
+ }