compressibleInterFoam: More consistent with interFoam and added partial support for LTS

applications/solvers/multiphase/compressibleInterFoam/Make/options

 EXE_INC = \
+    -I. \
+    -I../interFoam \
     -ItwoPhaseMixtureThermo \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \

applications/solvers/multiphase/compressibleInterFoam/TEqn.H

@@ -24,5 +24,6 @@
 
     fvOptions.correct(T);
 
+    mixture.correctThermo();
     mixture.correct();
 }

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -20,7 +20,7 @@
             fvc::reconstruct
             (
                 (
-                    interface.surfaceTensionForce()
+                    mixture.surfaceTensionForce()
                   - ghf*fvc::snGrad(rho)
                   - fvc::snGrad(p_rgh)
                 ) * mesh.magSf()

applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H → applications/solvers/multiphase/compressibleInterFoam/alphaEqn.H

@@ -2,48 +2,11 @@
     word alphaScheme("div(phi,alpha)");
     word alpharScheme("div(phirb,alpha)");
 
-    surfaceScalarField phir(phic*interface.nHatf());
+    surfaceScalarField phir(phic*mixture.nHatf());
 
     for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
     {
-        volScalarField::Internal Sp
-        (
-            IOobject
-            (
-                "Sp",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
-        );
-
-        volScalarField::Internal Su
-        (
-            IOobject
-            (
-                "Su",
-                runTime.timeName(),
-                mesh
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            divU*min(alpha1, scalar(1))
-        );
-
-        forAll(dgdt, celli)
-        {
-            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
-            {
-                Sp[celli] -= dgdt[celli]*alpha1[celli];
-                Su[celli] += dgdt[celli]*alpha1[celli];
-            }
-            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
-            {
-                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
-            }
-        }
-
+        #include "alphaSuSp.H"
 
         surfaceScalarField alphaPhi1
         (

applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H → applications/solvers/multiphase/compressibleInterFoam/alphaEqnSubCycle.H

 {
-    #include "alphaControls.H"
-
     surfaceScalarField phic(mag(phi/mesh.magSf()));
-    phic = min(interface.cAlpha()*phic, max(phic));
+    phic = min(mixture.cAlpha()*phic, max(phic));
 
     volScalarField divU(fvc::div(fvc::absolute(phi, U)));
 
@@ -27,7 +25,7 @@
             !(++alphaSubCycle).end();
         )
         {
-            #include "alphaEqns.H"
+            #include "alphaEqn.H"
             rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
         }
 
@@ -35,6 +33,6 @@
     }
     else
     {
-        #include "alphaEqns.H"
+        #include "alphaEqn.H"
     }
 }

applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H

+        volScalarField::Internal Sp
+        (
+            IOobject
+            (
+                "Sp",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
+        );
+
+        volScalarField::Internal Su
+        (
+            IOobject
+            (
+                "Su",
+                runTime.timeName(),
+                mesh
+            ),
+            // Divergence term is handled explicitly to be
+            // consistent with the explicit transport solution
+            divU*min(alpha1, scalar(1))
+        );
+
+        forAll(dgdt, celli)
+        {
+            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
+            {
+                Sp[celli] -= dgdt[celli]*alpha1[celli];
+                Su[celli] += dgdt[celli]*alpha1[celli];
+            }
+            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
+            {
+                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
+            }
+        }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options

 EXE_INC = \
     -I.. \
+    -I../../interFoam \
     -I../twoPhaseMixtureThermo \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \

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

@@ -41,13 +41,14 @@ Description
 #include "dynamicFvMesh.H"
 #include "MULES.H"
 #include "subCycle.H"
-#include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "twoPhaseMixtureThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
 #include "CorrectPhi.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,64 +82,67 @@ int main(int argc, char *argv[])
     {
         #include "readControls.H"
 
-        {
-            // Store divU from the previous mesh so that it can be mapped
-            // and used in correctPhi to ensure the corrected phi has the
-            // same divergence
-            volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
+        // Store divU from the previous mesh so that it can be mapped
+        // and used in correctPhi to ensure the corrected phi has the
+        // same divergence
+        volScalarField divU("divU0", fvc::div(fvc::absolute(phi, U)));
 
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
             #include "CourantNo.H"
+            #include "alphaCourantNo.H"
             #include "setDeltaT.H"
+        }
 
-            runTime++;
-
-            Info<< "Time = " << runTime.timeName() << nl << endl;
-
-            scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+        runTime++;
 
-            // Do any mesh changes
-            mesh.update();
+        Info<< "Time = " << runTime.timeName() << nl << endl;
 
-            if (mesh.changing())
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
             {
-                Info<< "Execution time for mesh.update() = "
-                    << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
-                    << " s" << endl;
+                scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-                gh = (g & mesh.C()) - ghRef;
-                ghf = (g & mesh.Cf()) - ghRef;
-            }
+                mesh.update();
 
-            if (mesh.changing() && correctPhi)
-            {
-                // Calculate absolute flux from the mapped surface velocity
-                phi = mesh.Sf() & Uf;
+                if (mesh.changing())
+                {
+                    Info<< "Execution time for mesh.update() = "
+                        << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                        << " s" << endl;
 
-                #include "correctPhi.H"
+                    gh = (g & mesh.C()) - ghRef;
+                    ghf = (g & mesh.Cf()) - ghRef;
+                }
 
-                // Make the fluxes relative to the mesh motion
-                fvc::makeRelative(phi, U);
-            }
-        }
+                if (mesh.changing() && correctPhi)
+                {
+                    // Calculate absolute flux from the mapped surface velocity
+                    phi = mesh.Sf() & Uf;
 
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
-        }
+                    #include "correctPhi.H"
 
-        turbulence->correct();
+                    // Make the fluxes relative to the mesh motion
+                    fvc::makeRelative(phi, U);
 
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            #include "alphaEqnsSubCycle.H"
+                    mixture.correct();
+                }
 
-            // correct interface on first PIMPLE corrector
-            if (pimple.corr() == 1)
-            {
-                interface.correct();
+                if (mesh.changing() && checkMeshCourantNo)
+                {
+                    #include "meshCourantNo.H"
+                }
             }
 
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             #include "UEqn.H"
@@ -149,13 +153,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
-
-        rho = alpha1*rho1 + alpha2*rho2;
 
-        // Correct p_rgh for consistency with p and the updated densities
-        p_rgh = p - rho*gh;
-        p_rgh.correctBoundaryConditions();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 

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

@@ -9,3 +9,8 @@ bool checkMeshCourantNo
 (
     pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false)
 );
+
+bool moveMeshOuterCorrectors
+(
+    pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false)
+);

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

@@ -12,7 +12,7 @@
     surfaceScalarField phig
     (
         (
-            interface.surfaceTensionForce()
+            mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf()
     );

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

@@ -4,3 +4,6 @@ correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
 
 checkMeshCourantNo =
     pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
+
+moveMeshOuterCorrectors =
+    pimple.dict().lookupOrDefault<Switch>("moveMeshOuterCorrectors", false);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -39,12 +39,13 @@ Description
 #include "MULES.H"
 #include "subCycle.H"
 #include "rhoThermo.H"
-#include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "twoPhaseMixtureThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,8 +60,6 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "createFields.H"
     #include "createFvOptions.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
 
     volScalarField& p = mixture.p();
     volScalarField& T = mixture.T();
@@ -69,6 +68,13 @@ int main(int argc, char *argv[])
 
     turbulence->validate();
 
+    if (!LTS)
+    {
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -76,8 +82,17 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "CourantNo.H"
-        #include "setDeltaT.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNo.H"
+            #include "alphaCourantNo.H"
+            #include "setDeltaT.H"
+        }
 
         runTime++;
 
@@ -86,13 +101,8 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            #include "alphaEqnsSubCycle.H"
-
-            // correct interface on first PIMPLE corrector
-            if (pimple.corr() == 1)
-            {
-                interface.correct();
-            }
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
 
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 

applications/solvers/multiphase/compressibleInterFoam/createFields.H

+#include "createRDeltaT.H"
+
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
 (
@@ -29,7 +31,7 @@ volVectorField U
 #include "createPhi.H"
 
 Info<< "Constructing twoPhaseMixtureThermo\n" << endl;
-twoPhaseMixtureThermo mixture(mesh);
+twoPhaseMixtureThermo mixture(U, phi);
 
 volScalarField& alpha1(mixture.alpha1());
 volScalarField& alpha2(mixture.alpha2());
@@ -89,9 +91,6 @@ volScalarField dgdt
     pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
 );
 
-// Construct interface from alpha1 distribution
-interfaceProperties interface(alpha1, U, mixture);
-
 // Construct compressible turbulence model
 autoPtr<compressible::turbulenceModel> turbulence
 (

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -12,7 +12,7 @@
     surfaceScalarField phig
     (
         (
-            interface.surfaceTensionForce()
+            mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
         )*rAUf*mesh.magSf()
     );

applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/Make/options

@@ -2,6 +2,7 @@ EXE_INC = \
     -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)/finiteVolume/lnInclude
 
 LIB_LIBS = \
@@ -9,4 +10,5 @@ LIB_LIBS = \
     -lfluidThermophysicalModels \
     -lspecie \
     -ltwoPhaseMixture \
+    -linterfaceProperties \
     -lfiniteVolume

applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.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
@@ -40,11 +40,13 @@ namespace Foam
 
 Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
 (
-    const fvMesh& mesh
+    const volVectorField& U,
+    const surfaceScalarField& phi
 )
 :
-    psiThermo(mesh, word::null),
-    twoPhaseMixture(mesh, *this),
+    psiThermo(U.mesh(), word::null),
+    twoPhaseMixture(U.mesh(), *this),
+    interfaceProperties(alpha1(), U, *this),
     thermo1_(nullptr),
     thermo2_(nullptr)
 {
@@ -58,8 +60,8 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
         T2.write();
     }
 
-    thermo1_ = rhoThermo::New(mesh, phase1Name());
-    thermo2_ = rhoThermo::New(mesh, phase2Name());
+    thermo1_ = rhoThermo::New(U.mesh(), phase1Name());
+    thermo2_ = rhoThermo::New(U.mesh(), phase2Name());
 
     thermo1_->validate(phase1Name(), "e");
     thermo2_->validate(phase2Name(), "e");
@@ -76,17 +78,23 @@ Foam::twoPhaseMixtureThermo::~twoPhaseMixtureThermo()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::twoPhaseMixtureThermo::correct()
+void Foam::twoPhaseMixtureThermo::correctThermo()
 {
     thermo1_->he() = thermo1_->he(p_, T_);
     thermo1_->correct();
 
     thermo2_->he() = thermo2_->he(p_, T_);
     thermo2_->correct();
+}
+
 
+void Foam::twoPhaseMixtureThermo::correct()
+{
     psi_ = alpha1()*thermo1_->psi() + alpha2()*thermo2_->psi();
     mu_ = alpha1()*thermo1_->mu() + alpha2()*thermo2_->mu();
     alpha_ = alpha1()*thermo1_->alpha() + alpha2()*thermo2_->alpha();
+
+    interfaceProperties::correct();
 }
 
 

applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,7 @@ SourceFiles
 #include "rhoThermo.H"
 #include "psiThermo.H"
 #include "twoPhaseMixture.H"
+#include "interfaceProperties.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,7 +53,8 @@ namespace Foam
 class twoPhaseMixtureThermo
 :
     public psiThermo,
-    public twoPhaseMixture
+    public twoPhaseMixture,
+    public interfaceProperties
 {
     // Private data
 
@@ -71,10 +73,11 @@ public:
 
     // Constructors
 
-        //- Construct from mesh
+        //- Construct from components
         twoPhaseMixtureThermo
         (
-            const fvMesh& mesh
+            const volVectorField& U,
+            const surfaceScalarField& phi
         );
 
 
@@ -104,7 +107,10 @@ public:
             return thermo2_();
         }
 
-        //- Update properties
+        //- Correct the thermodynamics of each phase
+        virtual void correctThermo();
+
+        //- Update mixture properties
         virtual void correct();
 
         //- Return true if the equation of state is incompressible