Merge branch 'integration-foundation' into 'develop'

Integration openfoam.org

See merge request !144

applications/solvers/combustion/fireFoam/createFieldRefs.H

@@ -1,4 +1,4 @@
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-filmModelType& surfaceFilm = tsurfaceFilm();
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();
 const label inertIndex(composition.species()[inertSpecie]);

applications/solvers/combustion/fireFoam/createSurfaceFilmModel.H

@@ -1,5 +1,6 @@
 Info<< "\nConstructing surface film model" << endl;
 
-typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
+autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
-
+(
-autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
+    regionModels::surfaceFilmModel::New(mesh, g)
+);

applications/solvers/electromagnetics/mhdFoam/createPhiB.H

@@ -9,7 +9,7 @@ IOobject phiBHeader
 
 surfaceScalarField* phiBPtr = nullptr;
 
-    if (phiBHeader.typeHeaderOk<surfaceScalarField>(true))
+if (phiBHeader.typeHeaderOk<surfaceScalarField>(true))
 {
     Info<< "Reading face flux ";
 

applications/solvers/heatTransfer/thermoFoam/setAlphaEff.H

@@ -11,7 +11,7 @@ IOobject turbulencePropertiesHeader
     false
 );
 
-if (turbulencePropertiesHeader.typeHeaderOk<IOdictionary>(false))
+if (turbulencePropertiesHeader.typeHeaderOk<IOdictionary>(true))
 {
     autoPtr<compressible::turbulenceModel> turbulence
     (

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C

@@ -77,6 +77,9 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        // Store the particle positions
+        kinematicCloud.storeGlobalPositions();
+
         mesh.update();
 
         // Calculate absolute flux from the mapped surface velocity

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/icoUncoupledKinematicParcelDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,8 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        kinematicCloud.storeGlobalPositions();
+
         mesh.update();
 
         U.correctBoundaryConditions();

applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files

@@ -1,3 +0,0 @@
-reactingParcelFilmFoam.C
-
-EXE = $(FOAM_APPBIN)/reactingParcelFilmFoam

applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H

@@ -1,34 +0,0 @@
-    MRF.correctBoundaryVelocity(U);
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U) + fvm::div(phi, U)
-      + MRF.DDt(rho, U)
-      + turbulence->divDevRhoReff(U)
-     ==
-        parcels.SU(U)
-      + fvOptions(rho, U)
-    );
-
-    UEqn.relax();
-
-    fvOptions.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-          ==
-            fvc::reconstruct
-            (
-                (
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                )*mesh.magSf()
-            )
-        );
-
-        fvOptions.correct(U);
-        K = 0.5*magSqr(U);
-    }

applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H

@@ -1,9 +0,0 @@
-Info<< "\nConstructing reacting cloud" << endl;
-basicReactingCloud parcels
-(
-    "reactingCloud1",
-    rho,
-    U,
-    g,
-    slgThermo
-);

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -1,141 +0,0 @@
-Info<< "Creating combustion model\n" << endl;
-
-autoPtr<combustionModels::psiCombustionModel> combustion
-(
-    combustionModels::psiCombustionModel::New(mesh)
-);
-
-psiReactionThermo& thermo = combustion->thermo();
-thermo.validate(args.executable(), "h", "e");
-
-SLGThermo slgThermo(mesh, thermo);
-
-basicSpecieMixture& composition = thermo.composition();
-PtrList<volScalarField>& Y = composition.Y();
-
-const word inertSpecie(thermo.lookup("inertSpecie"));
-if (!composition.species().found(inertSpecie))
-{
-    FatalIOErrorIn(args.executable().c_str(), thermo)
-        << "Inert specie " << inertSpecie << " not found in available species "
-        << composition.species()
-        << exit(FatalIOError);
-}
-
-Info<< "Creating field rho\n" << endl;
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
-
-volScalarField& p = thermo.p();
-
-Info<< "\nReading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "compressibleCreatePhi.H"
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-// Set the turbulence into the combustion model
-combustion->setTurbulence(turbulence());
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-// Force p_rgh to be consistent with p
-p_rgh = p - rho*gh;
-
-mesh.setFluxRequired(p_rgh.name());
-
-multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-forAll(Y, i)
-{
-    fields.add(Y[i]);
-}
-fields.add(thermo.he());
-
-IOdictionary additionalControlsDict
-(
-    IOobject
-    (
-        "additionalControls",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ_IF_MODIFIED,
-        IOobject::NO_WRITE
-    )
-);
-
-Switch solvePrimaryRegion
-(
-    additionalControlsDict.lookup("solvePrimaryRegion")
-);
-
-volScalarField Qdot
-(
-    IOobject
-    (
-        "Qdot",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
-);
-
-#include "createDpdt.H"
-
-#include "createK.H"
-
-#include "createMRF.H"
-#include "createClouds.H"
-#include "createRadiationModel.H"
-#include "createSurfaceFilmModel.H"

applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H

@@ -1,5 +0,0 @@
-Info<< "\nConstructing surface film model" << endl;
-
-typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
-
-autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -1,59 +0,0 @@
-rho = thermo.rho();
-
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-surfaceScalarField phiHbyA
-(
-    "phiHbyA",
-    (
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, phi)
-    )
-  + phig
-);
-
-MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-// Update the pressure BCs to ensure flux consistency
-constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-while (pimple.correctNonOrthogonal())
-{
-    fvScalarMatrix p_rghEqn
-    (
-        fvc::ddt(psi, rho)*gh
-      + fvc::div(phiHbyA)
-      + fvm::ddt(psi, p_rgh)
-      - fvm::laplacian(rhorAUf, p_rgh)
-     ==
-        parcels.Srho()
-      + surfaceFilm.Srho()
-      + fvOptions(psi, p_rgh, rho.name())
-    );
-
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-    if (pimple.finalNonOrthogonalIter())
-    {
-        phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
-        U.correctBoundaryConditions();
-        fvOptions.correct(U);
-    }
-}
-
-p = p_rgh + rho*gh;
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-K = 0.5*magSqr(U);
-
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}

applications/solvers/lagrangian/reactingParcelFoam/EEqn.H

@@ -19,6 +19,7 @@
      ==
         rho*(U&g)
       + parcels.Sh(he)
+      + surfaceFilm.Sh()
       + radiation->Sh(thermo, he)
       + Qdot
       + fvOptions(rho, he)
@@ -35,6 +36,6 @@
     thermo.correct();
     radiation->correct();
 
-    Info<< "T gas min/max   " << min(T).value() << ", "
+    Info<< "T gas min/max   = " << min(T).value() << ", "
         << max(T).value() << endl;
 }

applications/solvers/lagrangian/reactingParcelFoam/Make/options

@@ -1,12 +1,11 @@
 EXE_INC = \
     -I. \
+    -I../reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@@ -17,33 +16,33 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
-
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
+    -lsampling \
     -lmeshTools \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
-    -llagrangian \
-    -llagrangianIntermediate \
-    -llagrangianTurbulence \
     -lspecie \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lthermophysicalProperties \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lchemistryModel \
-    -lradiationModels \
-    -lODE \
     -lregionModels \
+    -lradiationModels \
     -lsurfaceFilmModels \
-    -lcombustionModels \
+    -lsurfaceFilmDerivedFvPatchFields \
-    -lfvOptions \
+    -llagrangian \
-    -lsampling
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lODE \
+    -lcombustionModels

applications/solvers/lagrangian/reactingParcelFoam/UEqn.H

@@ -6,8 +6,7 @@
       + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
-        rho()*g
+        parcels.SU(U)
-      + parcels.SU(U)
       + fvOptions(rho, U)
     );
 
@@ -17,7 +16,18 @@
 
     if (pimple.momentumPredictor())
     {
-        solve(UEqn == -fvc::grad(p));
+        solve
+        (
+            UEqn
+          ==
+            fvc::reconstruct
+            (
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            )
+        );
 
         fvOptions.correct(U);
         K = 0.5*magSqr(U);

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -9,6 +9,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     )
 );
 
+
 {
     combustion->correct();
     Qdot = combustion->Qdot();
@@ -24,11 +25,12 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
+              - fvm::laplacian(turbulence->alphaEff(), Yi)
-             ==
+              ==
                 parcels.SYi(i, Yi)
-              + combustion->R(Yi)
               + fvOptions(rho, Yi)
+              + combustion->R(Yi)
+              + surfaceFilm.Srho(i)
             );
 
             YEqn.relax();

applications/solvers/lagrangian/reactingParcelFoam/createFieldRefs.H

@@ -1,3 +1,5 @@
+const label inertIndex(composition.species()[inertSpecie]);
+
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
-const label inertIndex(composition.species()[inertSpecie]);
+regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm();

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -1,7 +1,5 @@
 #include "createRDeltaT.H"
 
-#include "readGravitationalAcceleration.H"
-
 Info<< "Creating combustion model\n" << endl;
 
 autoPtr<combustionModels::rhoCombustionModel> combustion
@@ -26,8 +24,7 @@ if (!composition.species().found(inertSpecie))
         << exit(FatalIOError);
 }
 
-volScalarField& p = thermo.p();
+Info<< "Creating field rho\n" << endl;
-
 volScalarField rho
 (
     IOobject
@@ -41,6 +38,8 @@ volScalarField rho
     thermo.rho()
 );
 
+volScalarField& p = thermo.p();
+
 Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
@@ -57,30 +56,6 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-mesh.setFluxRequired(p.name());
-
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -96,6 +71,31 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Force p_rgh to be consistent with p
+p_rgh = p - rho*gh;
+
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
+mesh.setFluxRequired(p_rgh.name());
+
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -105,6 +105,11 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
+Switch solvePrimaryRegion
+(
+    pimple.dict().lookupOrDefault<Switch>("solvePrimaryRegion", true)
+);
+
 volScalarField Qdot
 (
     IOobject
@@ -126,3 +131,4 @@ volScalarField Qdot
 #include "createMRF.H"
 #include "createRadiationModel.H"
 #include "createClouds.H"
+#include "createSurfaceFilmModel.H"

applications/solvers/lagrangian/reactingParcelFoam/createSurfaceFilmModel.H

@@ -0,0 +1,6 @@
+Info<< "\nConstructing surface film model" << endl;
+
+autoPtr<regionModels::surfaceFilmModel> tsurfaceFilm
+(
+    regionModels::surfaceFilmModel::New(mesh, g)
+);

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -1,4 +1,7 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
 // pressure solution
@@ -7,6 +10,9 @@ const volScalarField psip0(psi*p);
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
@@ -14,58 +20,68 @@ surfaceScalarField phiHbyA
         fvc::flux(rho*HbyA)
       + rhorAUf*fvc::ddtCorr(rho, U, phi)
     )
+  + phig
 );
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 // Update the pressure BCs to ensure flux consistency
-constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
-fvScalarMatrix pDDtEqn
+fvScalarMatrix p_rghDDtEqn
 (
-    fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
   + fvc::div(phiHbyA)
  ==
     parcels.Srho()
-  + fvOptions(psi, p, rho.name())
+  + surfaceFilm.Srho()
+  + fvOptions(psi, p_rgh, rho.name())
 );
 
 while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix pEqn
+    fvScalarMatrix p_rghEqn
     (
-        pDDtEqn
+        p_rghDDtEqn
-      - fvm::laplacian(rhorAUf, p)
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi = phiHbyA + pEqn.flux();
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
     }
 }
 
-p.relax();
+p = p_rgh + rho*gh;
 
 // Thermodynamic density update
 thermo.correctRho(psi*p - psip0);
 
-#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
+#include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U = HbyA - rAU*fvc::grad(p);
+if (pressureControl.limit(p))
-U.correctBoundaryConditions();
+{
-fvOptions.correct(U);
+    p.correctBoundaryConditions();
-K = 0.5*magSqr(U);
+    rho = thermo.rho();
+    p_rgh = p - rho*gh;
+}
+else if (pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
 
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
 }
-
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-
-Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,18 +29,20 @@ Group
 
 Description
     Transient solver for compressible, turbulent flow with a reacting,
-    multiphase particle cloud, and optional sources/constraints.
+    multiphase particle cloud, and surface film modelling.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
 #include "basicReactingMultiphaseCloud.H"
+#include "surfaceFilmModel.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
-#include "fvOptions.H"
 #include "SLGThermo.H"
+#include "fvOptions.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
@@ -76,10 +78,14 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
 
-        if (!LTS)
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
         {
             #include "compressibleCourantNo.H"
-            #include "setDeltaT.H"
+            #include "setMultiRegionDeltaT.H"
         }
 
         runTime++;
@@ -87,34 +93,36 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         parcels.evolve();
+        surfaceFilm.evolve();
 
-        if (LTS)
+        if (solvePrimaryRegion)
         {
-            #include "setRDeltaT.H"
+            if (pimple.nCorrPIMPLE() <= 1)
-        }
-
-        #include "rhoEqn.H"
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            #include "UEqn.H"
-            #include "YEqn.H"
-            #include "EEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
             {
-                #include "pEqn.H"
+                #include "rhoEqn.H"
             }
 
-            if (pimple.turbCorr())
+            // --- PIMPLE loop
+            while (pimple.loop())
             {
-                turbulence->correct();
+                #include "UEqn.H"
-            }
+                #include "YEqn.H"
-        }
+                #include "EEqn.H"
 
-        rho = thermo.rho();
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+
+            rho = thermo.rho();
+        }
 
         runTime.write();
 
@@ -123,7 +131,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info<< "End\n" << endl;
+    Info<< "End" << endl;
 
     return 0;
 }

applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,15 +36,13 @@ Description
       + fvc::div(phi)
       ==
         parcels.Srho(rho)
+      + surfaceFilm.Srho()
       + fvOptions(rho)
     );
 
     rhoEqn.solve();
 
     fvOptions.correct(rho);
-
-    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
-        << endl;
 }
 
 // ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,25 +33,21 @@ Description
 
 if (adjustTimeStep)
 {
-    if (CoNum == -GREAT)
+    const scalar maxDeltaTFact =
-    {
+        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
-        CoNum = SMALL;
+    const scalar deltaTFact =
-    }
+        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
-
-    const scalar TFactorFluid = maxCo/(CoNum + SMALL);
-    const scalar TFactorFilm = maxCo/(surfaceFilm.CourantNumber() + SMALL);
-
-    const scalar dt0 = runTime.deltaTValue();
 
     runTime.setDeltaT
     (
         min
         (
-            dt0*min(min(TFactorFluid, TFactorFilm), 1.2),
+            deltaTFact*runTime.deltaTValue(),
             maxDeltaT
         )
     );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
 }
 
-
 // ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/simpleReactingParcelFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,9 +24,6 @@ License
 Application
     simpleReactingParcelFoam
 
-Group
-    grpLagrangianSolvers
-
 Description
     Steady state solver for compressible, turbulent flow with reacting,
     multiphase particle clouds and optional sources/constraints.
@@ -38,6 +35,7 @@ Description
 #include "basicReactingMultiphaseCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
+#include "IOporosityModelList.H"
 #include "fvOptions.H"
 #include "SLGThermo.H"
 #include "simpleControl.H"

applications/solvers/lagrangian/sprayFoam/EEqn.H

@@ -19,7 +19,6 @@
      ==
         rho*(U&g)
       + parcels.Sh(he)
-      + surfaceFilm.Sh()
       + radiation->Sh(thermo, he)
       + Qdot
       + fvOptions(rho, he)
@@ -36,6 +35,6 @@
     thermo.correct();
     radiation->correct();
 
-    Info<< "T gas min/max   = " << min(T).value() << ", "
+    Info<< "T gas min/max   " << min(T).value() << ", "
         << max(T).value() << endl;
 }

applications/solvers/lagrangian/sprayFoam/YEqn.H

@@ -9,7 +9,6 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     )
 );
 
-
 {
     combustion->correct();
     Qdot = combustion->Qdot();
@@ -25,12 +24,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->alphaEff(), Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
-              ==
+             ==
                 parcels.SYi(i, Yi)
-              + fvOptions(rho, Yi)
               + combustion->R(Yi)
-              + surfaceFilm.Srho(i)
+              + fvOptions(rho, Yi)
             );
 
             YEqn.relax();

applications/solvers/lagrangian/sprayFoam/createFieldRefs.H

@@ -1,5 +1,3 @@
-const label inertIndex(composition.species()[inertSpecie]);
-
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
-filmModelType& surfaceFilm = tsurfaceFilm();
+const label inertIndex(composition.species()[inertSpecie]);

applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,9 @@ int main(int argc, char *argv[])
             // Store momentum to set rhoUf for introduced faces.
             volVectorField rhoU("rhoU", rho*U);
 
+            // Store the particle positions
+            parcels.storeGlobalPositions();
+
             // Do any mesh changes
             mesh.update();
 

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/files

@@ -0,0 +1,3 @@
+uncoupledKinematicParcelDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/uncoupledKinematicParcelDyMFoam

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options

@@ -1,46 +1,36 @@
 EXE_INC = \
-    -I. \
+    -I.. \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I${LIB_SRC}/sampling/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I${LIB_SRC}/meshTools/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
-    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 EXE_LIBS = \
-    -lfiniteVolume \
+    -llagrangian \
-    -lfvOptions \
-    -lsampling \
-    -lmeshTools \
-    -lturbulenceModels \
-    -lcompressibleTurbulenceModels \
-    -lspecie \
-    -lcompressibleTransportModels \
-    -lfluidThermophysicalModels \
-    -lreactionThermophysicalModels \
-    -lSLGThermo \
-    -lchemistryModel \
-    -lregionModels \
-    -lradiationModels \
-    -lsurfaceFilmModels \
-    -lsurfaceFilmDerivedFvPatchFields \
     -llagrangianIntermediate \
     -llagrangianTurbulence \
-    -lODE \
+    -lcompressibleTransportModels \
-    -lcombustionModels
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lradiationModels \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -ldynamicMesh \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/uncoupledKinematicParcelDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,91 +22,57 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    reactingParcelFilmFoam
+    uncoupledKinematicParcelDyMFoam
-
-Group
-    grpLagrangianSolvers
 
 Description
-    Transient solver for compressible, turbulent flow with a reacting,
+    Transient solver for the passive transport of a particle cloud.
-    multiphase particle cloud, and surface film modelling.
+
+    Uses a pre- calculated velocity field to evolve the cloud.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "basicReactingCloud.H"
+#include "basicKinematicCloud.H"
-#include "surfaceFilmModel.H"
-#include "psiCombustionModel.H"
-#include "radiationModel.H"
-#include "SLGThermo.H"
-#include "fvOptions.H"
-#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+    argList::addOption
+    (
+        "cloudName",
+        "name",
+        "specify alternative cloud name. default is 'kinematicCloud'"
+    );
+
+    #define NO_CONTROL
     #include "postProcess.H"
 
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createMesh.H"
+    #include "createDynamicFvMesh.H"
-    #include "createControl.H"
-    #include "createTimeControls.H"
     #include "createFields.H"
-    #include "createFieldRefs.H"
-    #include "createFvOptions.H"
     #include "compressibleCourantNo.H"
-    #include "setInitialDeltaT.H"
-    #include "initContinuityErrs.H"
-
-    turbulence->validate();
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (runTime.loop())
     {
-        #include "readTimeControls.H"
-        #include "compressibleCourantNo.H"
-        #include "setMultiRegionDeltaT.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        parcels.evolve();
+        kinematicCloud.storeGlobalPositions();
 
-        surfaceFilm.evolve();
+        mesh.update();
 
-        if (solvePrimaryRegion)
+        U.correctBoundaryConditions();
-        {
-            #include "rhoEqn.H"
 
-            // --- PIMPLE loop
+        Info<< "Evolving " << kinematicCloud.name() << endl;
-            while (pimple.loop())
+        kinematicCloud.evolve();
-            {
-                #include "UEqn.H"
-                #include "YEqn.H"
-                #include "EEqn.H"
-
-                // --- Pressure corrector loop
-                while (pimple.correct())
-                {
-                    #include "pEqn.H"
-                }
-
-                if (pimple.turbCorr())
-                {
-                    turbulence->correct();
-                }
-            }
-
-            rho = thermo.rho();
-        }
 
         runTime.write();
 
@@ -115,7 +81,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info<< "End" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H

@@ -118,27 +118,27 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
-        alphaPhi = talphaPhiUD();
+        alphaPhi10 = talphaPhi1UD();
 
-        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
+        if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
             MULES::correct
             (
                 alphac,
                 alpha1,
-                alphaPhi,
+                alphaPhi10,
-                talphaPhiCorr0.ref(),
+                talphaPhi1Corr0.ref(),
                 zeroField(), zeroField(),
                 1, 0
             );
 
-            alphaPhi += talphaPhiCorr0();
+            alphaPhi10 += talphaPhi1Corr0();
         }
 
         // Cache the upwind-flux
-        talphaPhiCorr0 = talphaPhiUD;
+        talphaPhi1Corr0 = talphaPhi1UD;
 
         alpha2 = 1.0 - alpha1;
 
@@ -152,7 +152,7 @@
 
         surfaceScalarField phir(phic*mixture.nHatf());
 
-        tmp<surfaceScalarField> talphaPhiUn
+        tmp<surfaceScalarField> talphaPhi1Un
         (
             fvc::flux
             (
@@ -170,15 +170,15 @@
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
             volScalarField alpha10("alpha10", alpha1);
 
             MULES::correct
             (
                 alphac,
                 alpha1,
-                talphaPhiUn(),
+                talphaPhi1Un(),
-                talphaPhiCorr.ref(),
+                talphaPhi1Corr.ref(),
                 Sp,
                 (-Sp*alpha1)(),
                 1,
@@ -188,24 +188,24 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                alphaPhi += talphaPhiCorr();
+                alphaPhi10 += talphaPhi1Corr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                alphaPhi += 0.5*talphaPhiCorr();
+                alphaPhi10 += 0.5*talphaPhi1Corr();
             }
         }
         else
         {
-            alphaPhi = talphaPhiUn;
+            alphaPhi10 = talphaPhi1Un;
 
             MULES::explicitSolve
             (
                 alphac,
                 alpha1,
                 phiCN,
-                alphaPhi,
+                alphaPhi10,
                 Sp,
                 (Su + divU*min(alpha1(), scalar(1)))(),
                 1,
@@ -220,12 +220,12 @@
 
     if (alphaApplyPrevCorr && MULESCorr)
     {
-        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
+        talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
-        talphaPhiCorr0.ref().rename("alphaPhiCorr0");
+        talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
     }
     else
     {
-        talphaPhiCorr0.clear();
+        talphaPhi1Corr0.clear();
     }
 
     if
@@ -235,19 +235,20 @@
     )
     {
         #include "rhofs.H"
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
     }
     else
     {
         if (ocCoeff > 0)
         {
             // Calculate the end-of-time-step alpha flux
-            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
+            alphaPhi10 =
+                (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
         }
 
         // Calculate the end-of-time-step mass flux
         #include "rhofs.H"
-        rhoPhi = alphaPhi*(rho1f - rho2f) + alphaPhic*rho2f;
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + alphaPhic*rho2f;
     }
 
     Info<< "Phase-1 volume fraction = "

applications/solvers/multiphase/VoF/alphaEqn.H

@@ -65,6 +65,14 @@
         phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
     }
 
+    // Add the optional shear compression contribution
+    if (scAlpha > 0)
+    {
+        phic +=
+            scAlpha*mag(mesh.delta() & fvc::interpolate(symm(fvc::grad(U))));
+    }
+
+
     surfaceScalarField::Boundary& phicBf =
         phic.boundaryFieldRef();
 
@@ -105,6 +113,8 @@
                 phiCN,
                 upwind<scalar>(mesh, phiCN)
             ).fvmDiv(phiCN, alpha1)
+       // - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
+       //           + fvc::div(phiCN), alpha1)
          ==
             Su + fvm::Sp(Sp + divU, alpha1)
         );
@@ -117,19 +127,19 @@
             << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
             << endl;
 
-        tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
+        tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
-        alphaPhi = talphaPhiUD();
+        alphaPhi10 = talphaPhi1UD();
 
-        if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
+        if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
         {
             Info<< "Applying the previous iteration compression flux" << endl;
-            MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);
+            MULES::correct(alpha1, alphaPhi10, talphaPhi1Corr0.ref(), 1, 0);
 
-            alphaPhi += talphaPhiCorr0();
+            alphaPhi10 += talphaPhi1Corr0();
         }
 
         // Cache the upwind-flux
-        talphaPhiCorr0 = talphaPhiUD;
+        talphaPhi1Corr0 = talphaPhi1UD;
 
         alpha2 = 1.0 - alpha1;
 
@@ -143,7 +153,7 @@
 
         surfaceScalarField phir(phic*mixture.nHatf());
 
-        tmp<surfaceScalarField> talphaPhiUn
+        tmp<surfaceScalarField> talphaPhi1Un
         (
             fvc::flux
             (
@@ -161,15 +171,15 @@
 
         if (MULESCorr)
         {
-            tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
+            tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
             volScalarField alpha10("alpha10", alpha1);
 
             MULES::correct
             (
                 geometricOneField(),
                 alpha1,
-                talphaPhiUn(),
+                talphaPhi1Un(),
-                talphaPhiCorr.ref(),
+                talphaPhi1Corr.ref(),
                 Sp,
                 (-Sp*alpha1)(),
                 1,
@@ -179,24 +189,24 @@
             // Under-relax the correction for all but the 1st corrector
             if (aCorr == 0)
             {
-                alphaPhi += talphaPhiCorr();
+                alphaPhi10 += talphaPhi1Corr();
             }
             else
             {
                 alpha1 = 0.5*alpha1 + 0.5*alpha10;
-                alphaPhi += 0.5*talphaPhiCorr();
+                alphaPhi10 += 0.5*talphaPhi1Corr();
             }
         }
         else
         {
-            alphaPhi = talphaPhiUn;
+            alphaPhi10 = talphaPhi1Un;
 
             MULES::explicitSolve
             (
                 geometricOneField(),
                 alpha1,
                 phiCN,
-                alphaPhi,
+                alphaPhi10,
                 Sp,
                 (Su + divU*min(alpha1(), scalar(1)))(),
                 1,
@@ -211,34 +221,37 @@
 
     if (alphaApplyPrevCorr && MULESCorr)
     {
-        talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
+        talphaPhi1Corr0 = alphaPhi10 - talphaPhi1Corr0;
-        talphaPhiCorr0.ref().rename("alphaPhiCorr0");
+        talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
     }
     else
     {
-        talphaPhiCorr0.clear();
+        talphaPhi1Corr0.clear();
     }
 
+    #include "rhofs.H"
+
     if
     (
         word(mesh.ddtScheme("ddt(rho,U)"))
      == fv::EulerDdtScheme<vector>::typeName
+     || word(mesh.ddtScheme("ddt(rho,U)"))
+     == fv::localEulerDdtScheme<vector>::typeName
     )
     {
-        #include "rhofs.H"
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phiCN*rho2f;
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phiCN*rho2f;
     }
     else
     {
         if (ocCoeff > 0)
         {
             // Calculate the end-of-time-step alpha flux
-            alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
+            alphaPhi10 =
+                (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
         }
 
         // Calculate the end-of-time-step mass flux
-        #include "rhofs.H"
+        rhoPhi = alphaPhi10*(rho1f - rho2f) + phi*rho2f;
-        rhoPhi = alphaPhi*(rho1f - rho2f) + phi*rho2f;
     }
 
     Info<< "Phase-1 volume fraction = "

applications/solvers/multiphase/VoF/createAlphaFluxes.H

@@ -1,20 +1,21 @@
-IOobject alphaPhiHeader
+IOobject alphaPhi10Header
 (
-    "alphaPhi",
+    "alphaPhi10",
     runTime.timeName(),
     mesh,
     IOobject::READ_IF_PRESENT,
     IOobject::AUTO_WRITE
 );
 
-const bool alphaRestart = alphaPhiHeader.typeHeaderOk<surfaceScalarField>(true);
+const bool alphaRestart =
+    alphaPhi10Header.typeHeaderOk<surfaceScalarField>(true);
 
 // MULES flux from previous time-step
-surfaceScalarField alphaPhi
+surfaceScalarField alphaPhi10
 (
-    alphaPhiHeader,
+    alphaPhi10Header,
     phi*fvc::interpolate(alpha1)
 );
 
 // MULES Correction
-tmp<surfaceScalarField> talphaPhiCorr0;
+tmp<surfaceScalarField> talphaPhi1Corr0;

applications/solvers/multiphase/VoF/setRDeltaT.H

@@ -73,8 +73,8 @@
         rDeltaT.ref() = max
         (
             rDeltaT(),
-            pos(alpha1Bar() - alphaSpreadMin)
+            pos0(alpha1Bar() - alphaSpreadMin)
-           *pos(alphaSpreadMax - alpha1Bar())
+           *pos0(alphaSpreadMax - alpha1Bar())
            *fvc::surfaceSum(mag(phi))()()
            /((2*maxAlphaCo)*mesh.V())
         );

applications/solvers/multiphase/compressibleInterFoam/Allwclean

@@ -5,5 +5,6 @@ wclean libso twoPhaseMixtureThermo
 wclean libso surfaceTensionModels
 wclean
 wclean compressibleInterDyMFoam
+wclean compressibleInterFilmFoam
 
 #------------------------------------------------------------------------------

applications/solvers/multiphase/compressibleInterFoam/Allwmake

@@ -9,5 +9,6 @@ wmake $targetType surfaceTensionModels
 
 wmake $targetType
 wmake $targetType compressibleInterDyMFoam
+wmake $targetType compressibleInterFilmFoam
 
 #------------------------------------------------------------------------------

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -1,8 +1,8 @@
 {
     fvScalarMatrix TEqn
     (
-        fvm::ddt(rho, T)
+        fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
-      + fvm::div(rhoPhi, T)
+      - fvm::Sp(contErr, T)
       - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
       + (
             divU*p

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -1,6 +1,7 @@
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      - fvm::Sp(contErr, U)
       + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==

applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H

@@ -24,14 +24,14 @@ volScalarField::Internal Su
 
 forAll(dgdt, celli)
 {
-    if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
+    if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]*alpha1[celli];
+        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]*alpha1[celli];
+        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
     }
-    else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
+    else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
+        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H

@@ -1,3 +1,69 @@
+tmp<surfaceScalarField> talphaPhi1(alphaPhi10);
+
+if (nAlphaSubCycles > 1)
 {
-    #include "alphaEqnSubCycle.H"
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+
+    talphaPhi1 = new surfaceScalarField
+    (
+        IOobject
+        (
+            "alphaPhi1",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", alphaPhi10.dimensions(), 0)
+    );
+
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", 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)*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)
+    )()
+);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -150,8 +150,6 @@ int main(int argc, char *argv[])
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 
-            solve(fvm::ddt(rho) + fvc::div(rhoPhi));
-
             #include "UEqn.H"
             #include "TEqn.H"
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -56,13 +56,29 @@
     }
     else
     {
+        #include "rhofs.H"
+
         p_rghEqnComp1 =
-            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+            pos(alpha1)
-          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
+           *(
+                (
+                    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 =
-            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+            pos(alpha2)
-          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
+           *(
+               (
+                   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
@@ -78,16 +94,21 @@
 
         solve
         (
-            (
+            p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
-                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
-              + (max(alpha2, scalar(0))/rho2)*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
@@ -109,15 +130,9 @@
 
     rho = alpha1*rho1 + alpha2*rho2;
 
-    p = max(p_rgh + rho*gh, pMin);
+    // Correct p_rgh for consistency with p and the updated densities
     p_rgh = p - rho*gh;
     p_rgh.correctBoundaryConditions();
 
-    dgdt =
-    (
-        pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
-      - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
-    );
-
     K = 0.5*magSqr(U);
 }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/files

@@ -0,0 +1,6 @@
+VoFPatchTransfer/VoFPatchTransfer.C
+VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C
+VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C
+compressibleInterFilmFoam.C
+
+EXE = $(FOAM_APPBIN)/compressibleInterFilmFoam

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options

@@ -0,0 +1,44 @@
+EXE_INC = \
+    -I. \
+    -I.. \
+    -I../../VoF \
+    -I../twoPhaseMixtureThermo \
+    -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)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
+
+EXE_LIBS = \
+    -ltwoPhaseMixtureThermo \
+    -ltwoPhaseSurfaceTension \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -ltwoPhaseMixture \
+    -ltwoPhaseProperties \
+    -linterfaceProperties \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lSLGThermo \
+    -lsurfaceFilmModels \
+    -lsurfaceFilmDerivedFvPatchFields \
+    -ldynamicMesh \
+    -lmeshTools \
+    -ldynamicFvMesh \
+    -lfiniteVolume \
+    -lfvOptions

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/TEqn.H

@@ -0,0 +1,36 @@
+{
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(rho, T) + fvm::div(rhoPhi, T)
+      - fvm::Sp(contErr, T)
+      - fvm::laplacian(mixture.alphaEff(turbulence->mut()), T)
+      + (
+            (
+                fvc::div(fvc::absolute(phi, U), p)
+              + fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
+            )()() - contErr*K
+        )
+       *(
+           alpha1/mixture.thermo1().Cv()
+         + alpha2/mixture.thermo2().Cv()
+        )()()
+     ==
+        fvOptions(rho, T)
+      + pos(Srho)
+       *(
+            surfaceFilm.Sh()()/mixture.thermo1().Cp()()
+          + surfaceFilm.Tref*Srho
+        )
+    );
+
+    TEqn.relax();
+
+    fvOptions.constrain(TEqn);
+
+    TEqn.solve();
+
+    fvOptions.correct(T);
+
+    mixture.correctThermo();
+    mixture.correct();
+}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C

@@ -0,0 +1,328 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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 "VoFPatchTransfer.H"
+#include "twoPhaseMixtureThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(VoFPatchTransfer, 0);
+addToRunTimeSelectionTable(transferModel, VoFPatchTransfer, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+VoFPatchTransfer::VoFPatchTransfer
+(
+    surfaceFilmRegionModel& film,
+    const dictionary& dict
+)
+:
+    transferModel(type(), film, dict),
+    deltaFactorToVoF_
+    (
+        coeffDict_.lookupOrDefault<scalar>("deltaFactorToVoF", 1.0)
+    ),
+    deltaFactorToFilm_
+    (
+        coeffDict_.lookupOrDefault<scalar>("deltaFactorToFilm", 0.5)
+    ),
+    alphaToVoF_
+    (
+        coeffDict_.lookupOrDefault<scalar>("alphaToVoF", 0.5)
+    ),
+    alphaToFilm_
+    (
+        coeffDict_.lookupOrDefault<scalar>("alphaToFilm", 0.1)
+    ),
+    transferRateCoeff_
+    (
+        coeffDict_.lookupOrDefault<scalar>("transferRateCoeff", 0.1)
+    )
+{
+    const polyBoundaryMesh& pbm = film.regionMesh().boundaryMesh();
+    patchIDs_.setSize
+    (
+        pbm.size() - film.regionMesh().globalData().processorPatches().size()
+    );
+
+    if (coeffDict_.found("patches"))
+    {
+        const wordReList patchNames(coeffDict_.lookup("patches"));
+        const labelHashSet patchSet = pbm.patchSet(patchNames);
+
+        Info<< "        applying to patches:" << nl;
+
+        label pidi = 0;
+        forAllConstIter(labelHashSet, patchSet, iter)
+        {
+            const label patchi = iter.key();
+            patchIDs_[pidi++] = patchi;
+            Info<< "            " << pbm[patchi].name() << endl;
+        }
+        patchIDs_.setSize(pidi);
+        patchTransferredMasses_.setSize(pidi, 0);
+    }
+    else
+    {
+        Info<< "            applying to all patches" << endl;
+
+        forAll(patchIDs_, patchi)
+        {
+            patchIDs_[patchi] = patchi;
+        }
+
+        patchTransferredMasses_.setSize(patchIDs_.size(), 0);
+    }
+
+    if (!patchIDs_.size())
+    {
+        FatalErrorInFunction
+            << "No patches selected"
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+VoFPatchTransfer::~VoFPatchTransfer()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void VoFPatchTransfer::correct
+(
+    scalarField& availableMass,
+    scalarField& massToTransfer
+)
+{
+    NotImplemented;
+}
+
+
+void VoFPatchTransfer::correct
+(
+    scalarField& availableMass,
+    scalarField& massToTransfer,
+    scalarField& energyToTransfer
+)
+{
+    // Do not correct if no patches selected
+    if (!patchIDs_.size()) return;
+
+    const scalarField& delta = film().delta();
+    const scalarField& rho = film().rho();
+    const scalarField& magSf = film().magSf();
+
+    const polyBoundaryMesh& pbm = film().regionMesh().boundaryMesh();
+
+
+    const twoPhaseMixtureThermo& thermo
+    (
+        film().primaryMesh().lookupObject<twoPhaseMixtureThermo>
+        (
+            twoPhaseMixtureThermo::dictName
+        )
+    );
+
+    const volScalarField& heVoF = thermo.thermo1().he();
+    const volScalarField& TVoF = thermo.thermo1().T();
+    const volScalarField CpVoF(thermo.thermo1().Cp());
+    const volScalarField& rhoVoF = thermo.thermo1().rho()();
+    const volScalarField& alphaVoF = thermo.alpha1();
+
+    forAll(patchIDs_, pidi)
+    {
+        const label patchi = patchIDs_[pidi];
+        label primaryPatchi = -1;
+
+        forAll(film().intCoupledPatchIDs(), i)
+        {
+            const label filmPatchi = film().intCoupledPatchIDs()[i];
+
+            if (filmPatchi == patchi)
+            {
+                primaryPatchi = film().primaryPatchIDs()[i];
+            }
+        }
+
+        if (primaryPatchi != -1)
+        {
+            scalarField deltaCoeffs
+            (
+                film().primaryMesh().boundary()[primaryPatchi].deltaCoeffs()
+            );
+            film().toRegion(patchi, deltaCoeffs);
+
+            scalarField hp(heVoF.boundaryField()[primaryPatchi]);
+            film().toRegion(patchi, hp);
+
+            scalarField Tp(TVoF.boundaryField()[primaryPatchi]);
+            film().toRegion(patchi, Tp);
+
+            scalarField Cpp(CpVoF.boundaryField()[primaryPatchi]);
+            film().toRegion(patchi, Cpp);
+
+            scalarField rhop(rhoVoF.boundaryField()[primaryPatchi]);
+            film().toRegion(patchi, rhop);
+
+            scalarField alphap(alphaVoF.boundaryField()[primaryPatchi]);
+            film().toRegion(patchi, alphap);
+
+            scalarField Vp
+            (
+                film().primaryMesh().boundary()[primaryPatchi]
+               .patchInternalField(film().primaryMesh().V())
+            );
+            film().toRegion(patchi, Vp);
+
+            const polyPatch& pp = pbm[patchi];
+            const labelList& faceCells = pp.faceCells();
+
+            // Accumulate the total mass removed from patch
+            scalar dMassPatch = 0;
+
+            forAll(faceCells, facei)
+            {
+                const label celli = faceCells[facei];
+
+                scalar dMass = 0;
+
+                if
+                (
+                    delta[celli] > 2*deltaFactorToVoF_/deltaCoeffs[facei]
+                 || alphap[facei] > alphaToVoF_
+                )
+                {
+                    dMass =
+                        transferRateCoeff_*delta[celli]*rho[celli]*magSf[celli];
+
+                    massToTransfer[celli] += dMass;
+                    energyToTransfer[celli] += dMass*film().hs()[celli];
+                }
+
+                if
+                (
+                    alphap[facei] > 0
+                 && delta[celli] > 0
+                 && delta[celli] < 2*deltaFactorToFilm_/deltaCoeffs[facei]
+                 && alphap[facei] < alphaToFilm_
+                )
+                {
+                    dMass =
+                        -transferRateCoeff_*alphap[facei]*rhop[facei]*Vp[facei];
+
+                    massToTransfer[celli] += dMass;
+                    energyToTransfer[celli] += dMass*hp[facei];
+                }
+
+                availableMass[celli] -= dMass;
+                dMassPatch += dMass;
+            }
+
+            patchTransferredMasses_[pidi] += dMassPatch;
+            addToTransferredMass(dMassPatch);
+        }
+    }
+
+    transferModel::correct();
+
+    if (writeTime())
+    {
+        scalarField patchTransferredMasses0
+        (
+            getModelProperty<scalarField>
+            (
+                "patchTransferredMasses",
+                scalarField(patchTransferredMasses_.size(), 0)
+            )
+        );
+
+        scalarField patchTransferredMassTotals(patchTransferredMasses_);
+        Pstream::listCombineGather
+        (
+            patchTransferredMassTotals,
+            plusEqOp<scalar>()
+        );
+        patchTransferredMasses0 += patchTransferredMassTotals;
+
+        setModelProperty<scalarField>
+        (
+            "patchTransferredMasses",
+            patchTransferredMasses0
+        );
+
+        patchTransferredMasses_ = 0;
+    }
+}
+
+
+void VoFPatchTransfer::patchTransferredMassTotals
+(
+    scalarField& patchMasses
+) const
+{
+    // Do not correct if no patches selected
+    if (!patchIDs_.size()) return;
+
+    scalarField patchTransferredMasses
+    (
+        getModelProperty<scalarField>
+        (
+            "patchTransferredMasses",
+            scalarField(patchTransferredMasses_.size(), 0)
+        )
+    );
+
+    scalarField patchTransferredMassTotals(patchTransferredMasses_);
+    Pstream::listCombineGather(patchTransferredMassTotals, plusEqOp<scalar>());
+
+    forAll(patchIDs_, pidi)
+    {
+        const label patchi = patchIDs_[pidi];
+        patchMasses[patchi] +=
+            patchTransferredMasses[pidi] + patchTransferredMassTotals[pidi];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.H

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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::regionModels::surfaceFilmModels::VoFPatchTransfer
+
+Description
+    Transfer mass between the film and the VoF in the continuous phase.
+
+SourceFiles
+    VoFPatchTransfer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef VoFPatchTransfer_H
+#define VoFPatchTransfer_H
+
+#include "transferModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class VoFPatchTransfer Declaration
+\*---------------------------------------------------------------------------*/
+
+class VoFPatchTransfer
+:
+    public transferModel
+{
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        VoFPatchTransfer(const VoFPatchTransfer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const VoFPatchTransfer&);
+
+
+protected:
+
+        //- Factor of the cell height above which the film is transferred
+        //  to the VoF
+        scalar deltaFactorToVoF_;
+
+        //- Factor of the cell height below which the VoF may be transferred
+        //  to the film
+        scalar deltaFactorToFilm_;
+
+        //- VoF limit above which all of the film is transferred to the VoF
+        scalar alphaToVoF_;
+
+        //- VoF limit below which the VoF may be transferred to the film
+        scalar alphaToFilm_;
+
+        //- Transfer rate coefficient
+        scalar transferRateCoeff_;
+
+        //- List of patch IDs at which the film is removed
+        labelList patchIDs_;
+
+        //- Transferred mass for each patch at which the film is removed
+        scalarField patchTransferredMasses_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("VoFPatchTransfer");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        VoFPatchTransfer(surfaceFilmRegionModel& film, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~VoFPatchTransfer();
+
+
+    // Member Functions
+
+        //- Correct
+        virtual void correct
+        (
+            scalarField& availableMass,
+            scalarField& massToTransfer
+        );
+
+        //- Correct kinematic and thermodynamic transfers
+        virtual void correct
+        (
+            scalarField& availableMass,
+            scalarField& massToTransfer,
+            scalarField& energyToTransfer
+        );
+
+        //- Accumulate the total mass injected for the patches into the
+        //  scalarField provided
+        virtual void patchTransferredMassTotals
+        (
+            scalarField& patchMasses
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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 "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
+)
+:
+    cellSetOption(sourceName, modelType, dict, mesh),
+    alphaSolidT_(Function1<scalar>::New("alphaSolidT", coeffs_)),
+    L_("L", dimEnergy/dimMass, coeffs_),
+    relax_(coeffs_.lookupOrDefault("relax", 0.9)),
+    Cu_(coeffs_.lookupOrDefault<scalar>("Cu", 100000)),
+    q_(coeffs_.lookupOrDefault("q", 0.001)),
+    alphaSolid_
+    (
+        IOobject
+        (
+            alphaSolidName(),
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("alpha1", dimless, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    curTimeIndex_(-1)
+{
+    fieldNames_.setSize(2);
+    fieldNames_[0] = "U";
+    fieldNames_[1] = "T";
+    applied_.setSize(fieldNames_.size(), false);
+}
+
+
+// * * * * * * * * * * * * * * * 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);
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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::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 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);
+
+        //- Disallow default bitwise copy construct
+        VoFSolidificationMeltingSource(const VoFSolidificationMeltingSource&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const VoFSolidificationMeltingSource&);
+
+
+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
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceIO.C

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

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceTemplates.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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 "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_));
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/compressibleInterFilmFoam.C

@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 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.
+
+    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 "rhoThermo.H"
+#include "twoPhaseMixture.H"
+#include "twoPhaseMixtureThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "SLGThermo.H"
+#include "surfaceFilmModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createFields.H"
+    #include "createAlphaFluxes.H"
+    #include "createSurfaceFilmModel.H"
+    #include "createFvOptions.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();
+
+    turbulence->validate();
+
+    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"
+
+            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();
+
+        Info<< "ExecutionTime = "
+            << runTime.elapsedCpuTime()
+            << " s\n\n" << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/createSurfaceFilmModel.H

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

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/pEqn.H

@@ -0,0 +1,129 @@
+{
+    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)
+      + 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())
+    {
+        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
+        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
+
+        p_rghEqnComp1 =
+            fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
+          + correction
+            (
+                psi1*fvm::ddt(p_rgh)
+              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
+        p_rghEqnComp1.ref().relax();
+
+        p_rghEqnComp2 =
+            fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+          + correction
+            (
+                psi2*fvm::ddt(p_rgh)
+              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+            );
+        deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
+        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);
+}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -108,9 +108,7 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "alphaControls.H"
-            #include "alphaEqnSubCycle.H"
+            #include "compressibleAlphaEqnSubCycle.H"
-
-            solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
             volScalarField divU(fvc::div(fvc::absolute(phi, U)));

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -87,10 +87,7 @@ surfaceScalarField rhoPhi
     fvc::interpolate(rho)*phi
 );
 
-volScalarField dgdt
+volScalarField dgdt(alpha1*fvc::div(phi));
-(
-    pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
-);
 
 // Construct compressible turbulence model
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -53,13 +53,29 @@
     }
     else
     {
+        #include "rhofs.H"
+
         p_rghEqnComp1 =
-            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+            pos(alpha1)
-          + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
+           *(
+                (
+                    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 =
-            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+            pos(alpha2)
-          + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
+           *(
+               (
+                   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
@@ -75,11 +91,7 @@
 
         solve
         (
-            (
+            p_rghEqnComp1() + p_rghEqnComp2() + p_rghEqnIncomp,
-                (max(alpha1, scalar(0))/rho1)*p_rghEqnComp1()
-              + (max(alpha2, scalar(0))/rho2)*p_rghEqnComp2()
-            )
-          + p_rghEqnIncomp,
             mesh.solver(p_rgh.select(pimple.finalInnerIter()))
         );
 
@@ -90,8 +102,8 @@
 
             dgdt =
             (
-                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
+                alpha1*(p_rghEqnComp2 & p_rgh)
-              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
+              - alpha2*(p_rghEqnComp1 & p_rgh)
             );
 
             phi = phiHbyA + p_rghEqnIncomp.flux();

applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C

@@ -63,8 +63,8 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
     thermo1_ = rhoThermo::New(U.mesh(), phase1Name());
     thermo2_ = rhoThermo::New(U.mesh(), phase2Name());
 
-    thermo1_->validate(phase1Name(), "e");
+    // thermo1_->validate(phase1Name(), "e");
-    thermo2_->validate(phase2Name(), "e");
+    // thermo2_->validate(phase2Name(), "e");
 
     correct();
 }

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -950,7 +950,7 @@ Foam::multiphaseMixtureThermo::nearInterface() const
     forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
     {
         tnearInt.ref() =
-            max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
+            max(tnearInt(), pos0(phase() - 0.01)*pos0(0.99 - phase()));
     }
 
     return tnearInt;

applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H

@@ -104,8 +104,10 @@
             )
             {
                 phase().dgdt() =
-                    pos(phase())
+                    pos0(phase())
                   *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho();
+
+                phasei++;
             }
 
             phi = phiHbyA + p_rghEqnIncomp.flux();

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -132,7 +132,7 @@ int main(int argc, char *argv[])
                     // if the mesh topology changed
                     if (mesh.topoChanging())
                     {
-                        talphaPhiCorr0.clear();
+                        talphaPhi1Corr0.clear();
                     }
 
                     gh = (g & mesh.C()) - ghRef;

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqn.H

@@ -183,8 +183,8 @@
         solve(fvm::ddt(alpha1) + fvc::div(alphaPhi1));
 
         // Create the diffusion coefficients for alpha2<->alpha3
-        volScalarField Dc23(D23*max(alpha3, scalar(0))*pos(alpha2));
+        volScalarField Dc23(D23*max(alpha3, scalar(0))*pos0(alpha2));
-        volScalarField Dc32(D23*max(alpha2, scalar(0))*pos(alpha3));
+        volScalarField Dc32(D23*max(alpha2, scalar(0))*pos0(alpha3));
 
         // Add the diffusive flux for alpha3->alpha2
         alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);

applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C

@@ -214,8 +214,8 @@ Foam::threePhaseInterfaceProperties::nearInterface() const
 {
     return max
     (
-        pos(mixture_.alpha1() - 0.01)*pos(0.99 - mixture_.alpha1()),
+        pos0(mixture_.alpha1() - 0.01)*pos0(0.99 - mixture_.alpha1()),
-        pos(mixture_.alpha2() - 0.01)*pos(0.99 - mixture_.alpha2())
+        pos0(mixture_.alpha2() - 0.01)*pos0(0.99 - mixture_.alpha2())
     );
 }
 

applications/solvers/multiphase/interFoam/overInterDyMFoam/interDyMFoam.C

@@ -141,7 +141,7 @@ int main(int argc, char *argv[])
                     // if the mesh topology changed
                     if (mesh.topoChanging())
                     {
-                        talphaPhiCorr0.clear();
+                        talphaPhi1Corr0.clear();
                     }
 
                     gh = (g & mesh.C()) - ghRef;

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C

@@ -87,7 +87,7 @@ Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const
     return Pair<tmp<volScalarField>>
     (
         mcCoeff_*sqr(limitedAlpha1)*(1.0 - limitedAlpha1)
-       *pos(p - pSat())/max(p - pSat(), 0.01*pSat()),
+       *pos0(p - pSat())/max(p - pSat(), 0.01*pSat()),
 
         (-mvCoeff_)*limitedAlpha1*neg(p - pSat())
     );

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C

@@ -83,7 +83,7 @@ Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotP() const
 
     return Pair<tmp<volScalarField>>
     (
-        mcCoeff_*(1.0 - limitedAlpha1)*pos(p - pSat()),
+        mcCoeff_*(1.0 - limitedAlpha1)*pos0(p - pSat()),
         (-mvCoeff_)*limitedAlpha1*neg(p - pSat())
     );
 }

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C

@@ -136,7 +136,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotP() const
 
     return Pair<tmp<volScalarField>>
     (
-        Cc_*(1.0 - limitedAlpha1)*pos(p - pSat())*apCoeff,
+        Cc_*(1.0 - limitedAlpha1)*pos0(p - pSat())*apCoeff,
 
         (-Cv_)*(1.0 + alphaNuc() - limitedAlpha1)*neg(p - pSat())*apCoeff
     );

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,13 +78,13 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
     volScalarField Cds
     (
         neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
-      + pos(Re - 1000)*0.44
+      + pos0(Re - 1000)*0.44
     );
 
     // Wen and Yu (1966)
     return
     (
-        pos(alpha2 - 0.8)
+        pos0(alpha2 - 0.8)
        *(0.75*Cds*phase2_.rho()*Ur*bp/d)
       + neg(alpha2 - 0.8)
        *(

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
     volScalarField Cds
     (
         neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
-      + pos(Re - 1000)*0.44
+      + pos0(Re - 1000)*0.44
     );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,7 +74,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
     volScalarField Cds
     (
         neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
-      + pos(Re - 1000)*0.44
+      + pos0(Re - 1000)*0.44
     );
 
     return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,7 +75,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
     volScalarField B
     (
         neg(alpha2 - 0.85)*(0.8*pow(alpha2, 1.28))
-      + pos(alpha2 - 0.85)*(pow(alpha2, 2.65))
+      + pos0(alpha2 - 0.85)*(pow(alpha2, 2.65))
     );
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
     volScalarField Cds
     (
         neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
-      + pos(Re - 1000)*0.44
+      + pos0(Re - 1000)*0.44
     );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -846,7 +846,8 @@ Foam::multiphaseSystem::nearInterface() const
 
     forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
     {
-        tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
+        tnearInt.ref() =
+            max(tnearInt(), pos0(iter() - 0.01)*pos0(0.99 - iter()));
     }
 
     return tnearInt;

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -254,8 +254,8 @@
             // dgdt =
 
             // (
-            //     pos(alpha2)*(pEqnComp2 & p)/rho2
+            //     pos0(alpha2)*(pEqnComp2 & p)/rho2
-            //   - pos(alpha1)*(pEqnComp1 & p)/rho1
+            //   - pos0(alpha1)*(pEqnComp1 & p)/rho1
             // );
 
             p_rgh.relax();

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -548,7 +548,8 @@ Foam::multiphaseMixture::nearInterface() const
 
     forAllConstIter(PtrDictionary<phase>, phases_, iter)
     {
-        tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
+        tnearInt.ref() =
+            max(tnearInt(), pos0(iter() - 0.01)*pos0(0.99 - iter()));
     }
 
     return tnearInt;

applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/files

@@ -17,6 +17,7 @@ saturationModels/Antoine/Antoine.C
 saturationModels/AntoineExtended/AntoineExtended.C
 saturationModels/ArdenBuck/ArdenBuck.C
 saturationModels/polynomial/polynomial.C
+saturationModels/function1/function1.C
 saturationModels/constantSaturationConditions/constantSaturationConditions.C
 
 LIB = $(FOAM_LIBBIN)/libreactingEulerianInterfacialCompositionModels

applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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 "function1.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace saturationModels
+{
+    defineTypeNameAndDebug(function1, 0);
+    addToRunTimeSelectionTable(saturationModel, function1, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::saturationModels::function1::function1(const dictionary& dict)
+:
+    saturationModel(),
+    function_
+    (
+        Function1<scalar>::New("function", dict)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::saturationModels::function1::~function1()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::saturationModels::function1::pSat
+(
+    const volScalarField& T
+) const
+{
+    NotImplemented;
+    return volScalarField::null();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::saturationModels::function1::pSatPrime
+(
+    const volScalarField& T
+) const
+{
+    NotImplemented;
+    return volScalarField::null();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::saturationModels::function1::lnPSat
+(
+    const volScalarField& T
+) const
+{
+    NotImplemented;
+    return volScalarField::null();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::saturationModels::function1::Tsat
+(
+    const volScalarField& p
+) const
+{
+    tmp<volScalarField> tTsat
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Tsat",
+                p.mesh().time().timeName(),
+                p.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            p.mesh(),
+            dimensionedScalar("zero", dimTemperature, 0)
+        )
+    );
+
+    volScalarField& Tsat = tTsat.ref();
+
+    forAll(Tsat, celli)
+    {
+        Tsat[celli] = function_->value(p[celli]);
+    }
+
+    volScalarField::Boundary& TsatBf = Tsat.boundaryFieldRef();
+
+    forAll(Tsat.boundaryField(), patchi)
+    {
+        scalarField& Tsatp = TsatBf[patchi];
+        const scalarField& pp = p.boundaryField()[patchi];
+
+        forAll(Tsatp, facei)
+        {
+            Tsatp[facei] = function_->value(pp[facei]);
+
+        }
+    }
+
+    return tTsat;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/function1/function1.H

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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::saturationModels::function1
+
+Description
+    Saturation vapour temperature in terms of
+    the vapour pressure (in Pa). The saturation temperature in Kelvins is
+    specified as a Foam::Function1 type, to enable use of, e.g. constant,
+    polynomial, table values.
+
+    Currently this class only provides \f$T_sat\f$, the inverse function to
+    return the vapour pressure for a given temperature are not implemented.
+
+    Examples:
+
+    \verbatim
+        type function1;
+        function polynomial
+        (
+            (308.0422       0)
+            (0.0015096      1)
+            (-1.61589e-8    2)
+            (1.114106e-13   3)
+            (-4.52216e-19   4)
+            (1.05192e-24    5)
+            (-1.2953e-30    6)
+            (6.5365e-37     7)
+        )
+    \endverbatim
+
+    \verbatim
+        type function1;
+
+        function csvFile;
+        functionCoeffs
+        {
+            nHeaderLine         1;
+            refColumn           0;
+            componentColumns    (1);
+            separator           ",";
+            mergeSeparators     no;
+            file                "filename.csv";
+            outOfBounds         clamp;
+            interpolationScheme linear;
+        };
+    \endverbatim
+
+SourceFiles
+    function1.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef function1_saturationModel_H
+#define function1_saturationModel_H
+
+#include "saturationModel.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace saturationModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class function1 Declaration
+\*---------------------------------------------------------------------------*/
+
+class function1
+:
+    public saturationModel
+{
+    // Private data
+
+        //- Saturation temperature as a function of pressure
+        autoPtr<Function1<scalar>> function_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("function1");
+
+    // Constructors
+
+        //- Construct from a dictionary
+        function1(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~function1();
+
+
+    // Member Functions
+
+        //- Saturation pressure
+        virtual tmp<volScalarField> pSat(const volScalarField& T) const;
+
+        //- Saturation pressure derivetive w.r.t. temperature
+        virtual tmp<volScalarField> pSatPrime(const volScalarField& T) const;
+
+        //- Natural log of the saturation pressure
+        virtual tmp<volScalarField> lnPSat(const volScalarField& T) const;
+
+        //- Saturation temperature
+        virtual tmp<volScalarField> Tsat(const volScalarField& p) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace saturationModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/aspectRatioModels/VakhrushevEfremov/VakhrushevEfremov.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,9 +71,9 @@ Foam::aspectRatioModels::VakhrushevEfremov::E() const
 
     return
         neg(Ta - scalar(1))*scalar(1)
-      + pos(Ta - scalar(1))*neg(Ta - scalar(39.8))
+      + pos0(Ta - scalar(1))*neg(Ta - scalar(39.8))
        *pow3(0.81 + 0.206*tanh(1.6 - 2*log10(max(Ta, scalar(1)))))
-      + pos(Ta - scalar(39.8))*0.24;
+      + pos0(Ta - scalar(39.8))*0.24;
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Beetstra/Beetstra.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -73,11 +73,12 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Beetstra::CdRe() const
         max(scalar(1) - pair_.dispersed(), pair_.continuous().residualAlpha())
     );
 
-    volScalarField Re(pair_.Re());
+    volScalarField Res(alpha2*pair_.Re());
-    volScalarField ReLim
+
+    volScalarField ResLim
     (
         "ReLim",
-        max(Re, residualRe_)
+        max(Res, residualRe_)
     );
 
     volScalarField F0
@@ -89,9 +90,9 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Beetstra::CdRe() const
     volScalarField F1
     (
         "F1",
-        0.413*Re/(24.0*sqr(alpha2))*(1.0/alpha2
+        0.413*Res/(24.0*sqr(alpha2))*(1.0/alpha2
-        + 3.0*alpha1*alpha2 + 8.4*pow(ReLim, -0.343))
+        + 3.0*alpha1*alpha2 + 8.4*pow(ResLim, -0.343))
-        /(1.0 + pow(10.0, 3*alpha1)*pow(ReLim, -(1.0 + 4.0*alpha1)/2.0))
+        /(1.0 + pow(10.0, 3*alpha1)*pow(ResLim, -(1.0 + 4.0*alpha1)/2.0))
     );
 
     return 24.0*alpha2*(F0 + F1);

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,7 +84,7 @@ Foam::tmp<Foam::volScalarField>
 Foam::dragModels::GidaspowErgunWenYu::CdRe() const
 {
     return
-        pos(pair_.continuous() - 0.8)*WenYu_->CdRe()
+        pos0(pair_.continuous() - 0.8)*WenYu_->CdRe()
       + neg(pair_.continuous() - 0.8)*Ergun_->CdRe();
 }
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,7 +74,7 @@ Foam::dragModels::GidaspowSchillerNaumann::CdRe() const
     volScalarField CdsRe
     (
         neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687))/alpha2
-      + pos(Re - 1000)*0.44*max(Re, residualRe_)
+      + pos0(Re - 1000)*0.44*max(Re, residualRe_)
     );
 
     return

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,7 +80,7 @@ Foam::dragModels::IshiiZuber::CdRe() const
     volScalarField ReM(Re*muc/muMix);
     volScalarField CdRe
     (
-        pos(1000 - ReM)*24.0*(scalar(1) + 0.15*pow(ReM, 0.687))
+        pos0(1000 - ReM)*24.0*(scalar(1) + 0.15*pow(ReM, 0.687))
       + neg(1000 - ReM)*0.44*ReM
     );
 
@@ -92,7 +92,7 @@ Foam::dragModels::IshiiZuber::CdRe() const
     volScalarField CdReEllipse(Ealpha*0.6666*sqrt(Eo)*Re);
 
     return
-        pos(CdReEllipse - CdRe)
+        pos0(CdReEllipse - CdRe)
        *min(CdReEllipse, Re*sqr(1 - pair_.dispersed())*2.66667)
       + neg(CdReEllipse - CdRe)*CdRe;
 }

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Lain/Lain.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,9 +66,9 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Lain::CdRe() const
 
     return
         neg(Re - 1.5)*16.0
-      + pos(Re - 1.5)*neg(Re - 80.0)*14.9*pow(Re, 0.22)
+      + pos0(Re - 1.5)*neg(Re - 80.0)*14.9*pow(Re, 0.22)
-      + pos(Re - 80.0)*neg(Re - 1500.0)*48*(1.0 - 2.21/sqrt(max(Re, SMALL)))
+      + pos0(Re - 80.0)*neg(Re - 1500.0)*48*(1.0 - 2.21/sqrt(max(Re, SMALL)))
-      + pos(Re - 1500.0)*2.61*Re;
+      + pos0(Re - 1500.0)*2.61*Re;
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,7 +67,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::CdRe() const
 
     return
         neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687))
-      + pos(Re - 1000)*0.44*max(Re, residualRe_);
+      + pos0(Re - 1000)*0.44*max(Re, residualRe_);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,7 +71,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::CdRe() const
     volScalarField B
     (
         neg(alpha2 - 0.85)*(0.8*pow(alpha2, 1.28))
-      + pos(alpha2 - 0.85)*(pow(alpha2, 2.65))
+      + pos0(alpha2 - 0.85)*(pow(alpha2, 2.65))
     );
     volScalarField Re(pair_.Re());
     volScalarField Vr

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,15 +50,6 @@ Foam::dragModels::Tenneti::Tenneti
 )
 :
     dragModel(dict, pair, registerObject),
-    SchillerNaumann_
-    (
-        new SchillerNaumann
-        (
-            dict,
-            pair,
-            false
-        )
-    ),
     residualRe_("residualRe", dimless, dict.lookup("residualRe"))
 {}
 
@@ -83,6 +74,14 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Tenneti::CdRe() const
         max(scalar(1) - pair_.dispersed(), pair_.continuous().residualAlpha())
     );
 
+    volScalarField Res(alpha2*pair_.Re());
+
+    volScalarField CdReIsolated
+    (
+        neg(Res - 1000)*24.0*(1.0 + 0.15*pow(Res, 0.687))
+      + pos0(Res - 1000)*0.44*max(Res, residualRe_)
+    );
+
     volScalarField F0
     (
         5.81*alpha1/pow3(alpha2) + 0.48*pow(alpha1, 1.0/3.0)/pow4(alpha2)
@@ -90,16 +89,14 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::Tenneti::CdRe() const
 
     volScalarField F1
     (
-        pow(alpha1, 3)*max(pair_.Re(), residualRe_)
+        pow3(alpha1)*Res*(0.95 + 0.61*pow3(alpha1)/sqr(alpha2))
-       *(0.95 + 0.61*pow3(alpha1)/sqr(alpha2))
     );
 
     // Tenneti et al. correlation includes the mean pressure drag.
     // This was removed here by multiplying F by alpha2 for consistency with
     // the formulation used in OpenFOAM
     return
-        SchillerNaumann_->CdRe()/(alpha2*max(pair_.Re(), residualRe_)) +
+        CdReIsolated + 24.0*sqr(alpha2)*(F0 + F1);
-        24.0*sqr(alpha2)*(F0 + F1);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,9 +69,6 @@ class Tenneti
 {
     // Private data
 
-        //- SchillerNaumann drag model
-        autoPtr<SchillerNaumann> SchillerNaumann_;
-
         //- Residual Reynolds Number
         const dimensionedScalar residualRe_;
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,7 +72,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::CdRe() const
     volScalarField CdsRes
     (
         neg(Res - 1000)*24.0*(1.0 + 0.15*pow(Res, 0.687))
-      + pos(Res - 1000)*0.44*max(Res, residualRe_)
+      + pos0(Res - 1000)*0.44*max(Res, residualRe_)
     );
 
     return

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/liftModels/TomiyamaLift/TomiyamaLift.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,8 +70,8 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::TomiyamaLift::Cl() const
 
     return
         neg(EoH - scalar(4))*min(0.288*tanh(0.121*pair_.Re()), f)
-      + pos(EoH - scalar(4))*neg(EoH - scalar(10.7))*f
+      + pos0(EoH - scalar(4))*neg(EoH - scalar(10.7))*f
-      + pos(EoH - scalar(10.7))*(-0.288);
+      + pos0(EoH - scalar(10.7))*(-0.288);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,7 +101,11 @@ Foam::turbulentDispersionModels::Burns::D() const
            *sqr(pair_.dispersed().d())
         )
        *pair_.continuous().rho()
-       *(1.0 + pair_.dispersed()/max(pair_.continuous(), residualAlpha_));
+       *pair_.dispersed()
+       *(
+           1.0/max(pair_.dispersed(), residualAlpha_)
+         + 1.0/max(pair_.continuous(), residualAlpha_)
+        );
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,9 +80,9 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Frank::Fi() const
     return zeroGradWalls
     (
         (
-            pos(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179)
+            pos0(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179)
-          + pos(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187)
+          + pos0(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187)
-          + pos(Eo - 33.0)*0.179
+          + pos0(Eo - 33.0)*0.179
         )
        *max
         (

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,9 +78,9 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::Fi() const
     return zeroGradWalls
     (
         (
-            pos(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179)
+            pos0(Eo - 1.0)*neg(Eo - 5.0)*exp(-0.933*Eo + 0.179)
-          + pos(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187)
+          + pos0(Eo - 5.0)*neg(Eo - 33.0)*(0.00599*Eo - 0.0187)
-          + pos(Eo - 33.0)*0.179
+          + pos0(Eo - 33.0)*0.179
         )
        *0.5
        *pair_.dispersed().d()

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -379,8 +379,8 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
         (
             min
             (
-                (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
+                (pos0(iDmdt)*he2 + neg(iDmdt)*hef2)
-              - (neg(iDmdt)*he1 + pos(iDmdt)*hef1),
+              - (neg(iDmdt)*he1 + pos0(iDmdt)*hef1),
                 0.3*mag(hef2 - hef1)
             )
         );

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,6 +60,13 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::~AnisothermalPhaseModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class BasePhaseModel>
+bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
+{
+    return !this->thermo().incompressible();
+}
+
+
 template<class BasePhaseModel>
 void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctKinematics()
 {
@@ -135,7 +142,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
             he
         )
      ==
-        this->Qdot()
+        alpha*this->Qdot()
     );
 
     // Add the appropriate pressure-work term
@@ -156,13 +163,6 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 }
 
 
-template<class BasePhaseModel>
-bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
-{
-    return !this->thermo().incompressible();
-}
-
-
 template<class BasePhaseModel>
 const Foam::volScalarField&
 Foam::AnisothermalPhaseModel<BasePhaseModel>::K() const

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,6 +85,9 @@ public:
 
     // Member Functions
 
+        //- Return true if the phase is compressible otherwise false
+        virtual bool compressible() const;
+
         //- Correct the kinematics
         virtual void correctKinematics();
 
@@ -94,14 +97,8 @@ public:
         //- Return the enthalpy equation
         virtual tmp<fvScalarMatrix> heEqn();
 
-
+        //- Return the phase kinetic energy
-        // Compressibility (variable density)
+        virtual const volScalarField& K() const;
-
-            //- Return true if the phase is compressible otherwise false
-            virtual bool compressible() const;
-
-            //- Return the phase kinetic energy
-            virtual const volScalarField& K() const;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,6 +49,13 @@ Foam::IsothermalPhaseModel<BasePhaseModel>::~IsothermalPhaseModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class BasePhaseModel>
+bool Foam::IsothermalPhaseModel<BasePhaseModel>::compressible() const
+{
+    return !this->thermo().incompressible();
+}
+
+
 template<class BasePhaseModel>
 void Foam::IsothermalPhaseModel<BasePhaseModel>::correctThermo()
 {}

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,6 +71,9 @@ public:
 
     // Member Functions
 
+        //- Return true if the phase is compressible otherwise false
+        virtual bool compressible() const;
+
         //- Correct the thermodynamics
         virtual void correctThermo();
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -172,7 +172,7 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::YiEqn
             Yi
         )
      ==
-        this->R(Yi)
+        alpha*this->R(Yi)
 
       + fvc::ddt(residualAlpha_*rho, Yi)
       - fvm::ddt(residualAlpha_*rho, Yi)