Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

applications/solvers/combustion/chemFoam/createFields.H

@@ -12,7 +12,7 @@
         (
             "initialConditions",
             runTime.constant(),
-            runTime,
+            mesh,
             IOobject::MUST_READ_IF_MODIFIED,
             IOobject::NO_WRITE
         )
@@ -37,7 +37,7 @@
         (
             "rho",
             runTime.timeName(),
-            runTime,
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -50,7 +50,7 @@
         (
             "Rspecific",
             runTime.timeName(),
-            runTime,
+            mesh,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
@@ -69,7 +69,7 @@
         (
             "U",
             runTime.timeName(),
-            runTime,
+            mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),

applications/solvers/compressible/rhoCentralFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso BCs
 wclean

applications/solvers/compressible/rhoCentralFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 (wmake $targetType BCs && wmake $targetType && wmake $targetType rhoCentralDyMFoam)
 

applications/solvers/lagrangian/DPMFoam/Allwclean

@@ -1,7 +1,6 @@
 #!/bin/sh
 
 cd ${0%/*} || exit 1
-set -x
 
 wclean libso DPMTurbulenceModels
 wclean

applications/solvers/lagrangian/DPMFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType DPMTurbulenceModels
 

applications/solvers/multiphase/compressibleInterFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso twoPhaseMixtureThermo
 wclean

applications/solvers/multiphase/compressibleInterFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType twoPhaseMixtureThermo
 

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso multiphaseMixtureThermo
 wclean

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType multiphaseMixtureThermo
 wmake $targetType

applications/solvers/multiphase/driftFluxFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso mixtureViscosityModels
 wclean libso relativeVelocityModels

applications/solvers/multiphase/driftFluxFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType mixtureViscosityModels
 wmake $targetType relativeVelocityModels

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C

@@ -24,6 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "relativeVelocityModel.H"
+#include "fixedValueFvPatchFields.H"
+#include "slipFvPatchFields.H"
+#include "partialSlipFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -33,6 +36,34 @@ namespace Foam
     defineRunTimeSelectionTable(relativeVelocityModel, dictionary);
 }
 
+// * * * * * * * * * * * * * Private Member Functions   * * * * * * * * * * * //
+
+Foam::wordList Foam::relativeVelocityModel::UdmPatchFieldTypes() const
+{
+    const volVectorField& U = mixture_.U();
+
+    wordList UdmTypes
+    (
+        U.boundaryField().size(),
+        calculatedFvPatchScalarField::typeName
+    );
+
+    forAll(U.boundaryField(), i)
+    {
+        if
+        (
+            isA<fixedValueFvPatchVectorField>(U.boundaryField()[i])
+         || isA<slipFvPatchVectorField>(U.boundaryField()[i])
+         || isA<partialSlipFvPatchVectorField>(U.boundaryField()[i])
+        )
+        {
+            UdmTypes[i] = fixedValueFvPatchVectorField::typeName;
+        }
+    }
+
+    return UdmTypes;
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -55,11 +86,12 @@ Foam::relativeVelocityModel::relativeVelocityModel
             "Udm",
             alphac_.time().timeName(),
             alphac_.mesh(),
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         alphac_.mesh(),
-        dimensionedVector("Udm", dimVelocity, Zero)
+        dimensionedVector("Udm", dimVelocity, Zero),
+        UdmPatchFieldTypes()
     )
 {}
 

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,12 +52,16 @@ class relativeVelocityModel
 {
     // Private Member Functions
 
+        //- Return the list of patchFieldTypes for Udm derived from U
+        wordList UdmPatchFieldTypes() const;
+
         //- Disallow default bitwise copy construct
         relativeVelocityModel(const relativeVelocityModel&);
 
         //- Disallow default bitwise assignment
         void operator=(const relativeVelocityModel&);
 
+
 protected:
 
     // Protected data

applications/solvers/multiphase/interCondensingEvaporatingFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso temperaturePhaseChangeTwoPhaseMixtures
 wclean

applications/solvers/multiphase/interPhaseChangeFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso phaseChangeTwoPhaseMixtures
 wclean

applications/solvers/multiphase/interPhaseChangeFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType phaseChangeTwoPhaseMixtures
 wmake $targetType

applications/solvers/multiphase/multiphaseEulerFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso multiphaseSystem
 wclean libso interfacialModels

applications/solvers/multiphase/multiphaseEulerFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmakeLnInclude interfacialModels
 wmake $targetType multiphaseSystem

applications/solvers/multiphase/multiphaseInterFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso multiphaseMixture
 wclean

applications/solvers/multiphase/multiphaseInterFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType multiphaseMixture
 wmake $targetType

applications/solvers/multiphase/reactingEulerFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso phaseSystems
 wclean libso interfacialModels

applications/solvers/multiphase/reactingEulerFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmakeLnInclude interfacialModels
 wmakeLnInclude interfacialCompositionModels

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

@@ -77,11 +77,39 @@ void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctThermo()
 }
 
 
+template<class BasePhaseModel>
+Foam::tmp<Foam::volScalarField>
+Foam::AnisothermalPhaseModel<BasePhaseModel>::filterPressureWork
+(
+    const tmp<volScalarField>& pressureWork
+) const
+{
+    const volScalarField& alpha = *this;
+
+    scalar pressureWorkAlphaLimit =
+        this->thermo_->lookupOrDefault("pressureWorkAlphaLimit", 0.0);
+
+    if (pressureWorkAlphaLimit > 0)
+    {
+        return
+        (
+            max(alpha - pressureWorkAlphaLimit, scalar(0))
+           /max(alpha - pressureWorkAlphaLimit, pressureWorkAlphaLimit)
+        )*pressureWork;
+    }
+    else
+    {
+        return pressureWork;
+    }
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 {
     const volScalarField& alpha = *this;
+    const volVectorField& U = this->U();
     const surfaceScalarField& alphaPhi = this->alphaPhi();
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi();
 
@@ -93,7 +121,8 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 
     tmp<fvScalarMatrix> tEEqn
     (
-        fvm::ddt(alpha, this->rho(), he) + fvm::div(alphaRhoPhi, he)
+        fvm::ddt(alpha, this->rho(), he)
+      + fvm::div(alphaRhoPhi, he)
       - fvm::Sp(contErr, he)
 
       + fvc::ddt(alpha, this->rho(), K_) + fvc::div(alphaRhoPhi, K_)
@@ -112,13 +141,15 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
     // Add the appropriate pressure-work term
     if (he.name() == this->thermo_->phasePropertyName("e"))
     {
-        tEEqn.ref() +=
-            fvc::ddt(alpha)*this->thermo().p()
-          + fvc::div(alphaPhi, this->thermo().p());
+        tEEqn.ref() += filterPressureWork
+        (
+            fvc::div(fvc::absolute(alphaPhi, alpha, U), this->thermo().p())
+          + this->thermo().p()*fvc::ddt(alpha)
+        );
     }
     else if (this->thermo_->dpdt())
     {
-        tEEqn.ref() -= alpha*this->fluid().dpdt();
+        tEEqn.ref() -= filterPressureWork(alpha*this->fluid().dpdt());
     }
 
     return tEEqn;

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

@@ -58,6 +58,15 @@ class AnisothermalPhaseModel
         volScalarField K_;
 
 
+    // Private member functions
+
+        //- Optionally filter the pressure work term as the phase-fraction -> 0
+        tmp<volScalarField> filterPressureWork
+        (
+            const tmp<volScalarField>& pressureWork
+        ) const;
+
+
 public:
 
     // Constructors

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ReactingPhaseModel/ReactingPhaseModel.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-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,12 @@ Foam::ReactingPhaseModel<BasePhaseModel, ReactionType>::ReactingPhaseModel
     BasePhaseModel(fluid, phaseName, index, false),
     reaction_
     (
-        ReactionType::New(fluid.mesh(), this->name())
+        ReactionType::New
+        (
+            fluid.mesh(),
+            combustionModel::combustionPropertiesName,
+            this->name()
+        )
     )
 {
     this->thermo_ = &reaction_->thermo();

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso multiphaseSystem
 wclean libso multiphaseCompressibleTurbulenceModels

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType multiphaseSystem
 wmake $targetType multiphaseCompressibleTurbulenceModels

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H

@@ -35,6 +35,7 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
     }
 
     fluid.correctThermo();
+    fluid.correct();
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H

@@ -89,9 +89,6 @@ PtrList<surfaceScalarField> phiFs(phases.size());
 // --- Pressure corrector loop
 while (pimple.correct())
 {
-    // Update continuity errors due to temperature changes
-    fluid.correct();
-
     volScalarField rho("rho", fluid.rho());
 
     // Correct p_rgh for consistency with p and the updated densities

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C

@@ -77,11 +77,6 @@ int main(int argc, char *argv[])
 
     //#include "pUf/createDDtU.H"
 
-    int nEnergyCorrectors
-    (
-        pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1)
-    );
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -90,6 +85,11 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
 
+        int nEnergyCorrectors
+        (
+            pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1)
+        );
+
         if (LTS)
         {
             #include "setRDeltaT.H"

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso twoPhaseSystem
 wclean libso twoPhaseCompressibleTurbulenceModels

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmake $targetType twoPhaseSystem
 wmake $targetType twoPhaseCompressibleTurbulenceModels

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/EEqns.H

@@ -29,36 +29,25 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
     }
 
     {
-        tmp<fvScalarMatrix> E2eqn(phase2.heEqn());
+        tmp<fvScalarMatrix> E2Eqn(phase2.heEqn());
 
-        if (E2eqn.valid())
+        if (E2Eqn.valid())
         {
-            E2eqn =
+            E2Eqn =
             (
-                E2eqn
+                E2Eqn
              ==
                *heatTransfer[phase2.name()]
               + alpha2*rho2*(U2&g)
               + fvOptions(alpha2, rho2, phase2.thermo().he())
             );
 
-            E2eqn->relax();
-            fvOptions.constrain(E2eqn.ref());
-            E2eqn->solve();
+            E2Eqn->relax();
+            fvOptions.constrain(E2Eqn.ref());
+            E2Eqn->solve();
         }
     }
 
     fluid.correctThermo();
+    fluid.correct();
 }
-
-Info<< phase1.name() << " min/max T "
-    << min(phase1.thermo().T()).value()
-    << " - "
-    << max(phase1.thermo().T()).value()
-    << endl;
-
-Info<< phase2.name() << " min/max T "
-    << min(phase2.thermo().T()).value()
-    << " - "
-    << max(phase2.thermo().T()).value()
-    << endl;

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createRDeltaTf.H

@@ -0,0 +1,21 @@
+tmp<surfaceScalarField> trDeltaTf;
+
+if (LTS && faceMomentum)
+{
+    trDeltaTf = tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                fv::localEulerDdt::rDeltaTfName,
+                runTime.timeName(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("one", dimless/dimTime, 1)
+        )
+    );
+}

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -83,9 +83,6 @@ tmp<surfaceScalarField> phiF2;
 // --- Pressure corrector loop
 while (pimple.correct())
 {
-    // Update continuity errors due to temperature changes
-    fluid.correct();
-
     volScalarField rho("rho", fluid.rho());
 
     // Correct p_rgh for consistency with p and the updated densities

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H

@@ -1,9 +1,2 @@
-    ddtPhi1 =
-    (
-        (phi1 - phi1.oldTime())/runTime.deltaT()
-    );
-
-    ddtPhi2 =
-    (
-        (phi2 - phi2.oldTime())/runTime.deltaT()
-    );
+ddtPhi1 = fvc::ddt(phi1);
+ddtPhi2 = fvc::ddt(phi2);

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H

@@ -1,9 +1,2 @@
-    surfaceScalarField ddtPhi1
-    (
-        (phi1 - phi1.oldTime())/runTime.deltaT()
-    );
-
-    surfaceScalarField ddtPhi2
-    (
-        (phi2 - phi2.oldTime())/runTime.deltaT()
-    );
+surfaceScalarField ddtPhi1(fvc::ddt(phi1));
+surfaceScalarField ddtPhi2(fvc::ddt(phi2));

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -32,10 +32,10 @@ surfaceScalarField rAUf1
     IOobject::groupName("rAUf", phase1.name()),
     1.0
    /(
-        (alphaRhof10 + Vmf)/runTime.deltaT()
+        byDt(alphaRhof10 + Vmf)
       + fvc::interpolate(U1Eqn.A())
       + Kdf
-   )
+    )
 );
 
 surfaceScalarField rAUf2
@@ -43,10 +43,10 @@ surfaceScalarField rAUf2
     IOobject::groupName("rAUf", phase2.name()),
     1.0
    /(
-        (alphaRhof20 + Vmf)/runTime.deltaT()
+        byDt(alphaRhof20 + Vmf)
       + fvc::interpolate(U2Eqn.A())
       + Kdf
-   )
+    )
 );
 
 
@@ -93,9 +93,6 @@ tmp<surfaceScalarField> Ff2;
 
 while (pimple.correct())
 {
-    // Update continuity errors due to temperature changes
-    fluid.correct();
-
     volScalarField rho("rho", fluid.rho());
 
     // Correct p_rgh for consistency with p and the updated densities
@@ -160,7 +157,7 @@ while (pimple.correct())
         rAUf1
        *(
             (alphaRhof10 + Vmf)
-           *MRF.absolute(phi1.oldTime())/runTime.deltaT()
+           *byDt(MRF.absolute(phi1.oldTime()))
           + fvc::flux(U1Eqn.H())
           + Vmf*ddtPhi2
           + Kdf*MRF.absolute(phi2)
@@ -178,7 +175,7 @@ while (pimple.correct())
         rAUf2
        *(
             (alphaRhof20 + Vmf)
-           *MRF.absolute(phi2.oldTime())/runTime.deltaT()
+           *byDt(MRF.absolute(phi2.oldTime()))
           + fvc::flux(U2Eqn.H())
           + Vmf*ddtPhi1
           + Kdf*MRF.absolute(phi1)

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C

@@ -43,6 +43,22 @@ Description
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
+namespace Foam
+{
+    tmp<surfaceScalarField> byDt(const surfaceScalarField& sf)
+    {
+        if (fv::localEulerDdt::enabled(sf.mesh()))
+        {
+            return fv::localEulerDdt::localRDeltaTf(sf.mesh())*sf;
+        }
+        else
+        {
+            return sf/sf.mesh().time().deltaT();
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -76,13 +92,9 @@ int main(int argc, char *argv[])
         )
     );
 
+    #include "createRDeltaTf.H"
     #include "pUf/createDDtU.H"
 
-    int nEnergyCorrectors
-    (
-        pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1)
-    );
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -91,9 +103,18 @@ int main(int argc, char *argv[])
     {
         #include "readTimeControls.H"
 
+        int nEnergyCorrectors
+        (
+            pimple.dict().lookupOrDefault<int>("nEnergyCorrectors", 1)
+        );
+
         if (LTS)
         {
             #include "setRDeltaT.H"
+            if (faceMomentum)
+            {
+                #include "setRDeltaTf.H"
+            }
         }
         else
         {

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/setRDeltaTf.H

@@ -0,0 +1 @@
+trDeltaTf.ref() = fvc::interpolate(fv::localEulerDdt::localRDeltaT(mesh));

applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean

@@ -1,6 +1,5 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
-set -x
 
 wclean libso twoPhaseSystem
 wclean libso interfacialModels

applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-set -x
 
 wmakeLnInclude interfacialModels
 wmake $targetType twoPhaseSystem

applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H

@@ -16,7 +16,8 @@
       - contErr1*K1
       + (
             he1.name() == thermo1.phasePropertyName("e")
-          ? fvc::ddt(alpha1)*p + fvc::div(alphaPhi1, p)
+          ? fvc::div(fvc::absolute(alphaPhi1, alpha1, U1), p)
+          + p*fvc::ddt(alpha1)
           : -alpha1*dpdt
         )
 
@@ -48,7 +49,8 @@
       - contErr2*K2
       + (
             he2.name() == thermo2.phasePropertyName("e")
-          ? fvc::ddt(alpha2)*p + fvc::div(alphaPhi2, p)
+          ? fvc::div(fvc::absolute(alphaPhi2, alpha2, U2), p)
+          + p*fvc::ddt(alpha1)
           : -alpha2*dpdt
         )
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C

@@ -248,4 +248,30 @@ bool Foam::phaseModel::read(const dictionary& phaseProperties)
 }
 
 
+void Foam::phaseModel::correctInflowFlux(surfaceScalarField& alphaPhi) const
+{
+    surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+
+    // Ensure that the flux at inflow BCs is preserved
+    forAll(alphaPhiBf, patchi)
+    {
+        fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
+
+        if (!alphaPhip.coupled())
+        {
+            const scalarField& phip = phi().boundaryField()[patchi];
+            const scalarField& alphap = boundaryField()[patchi];
+
+            forAll(alphaPhip, facei)
+            {
+                if (phip[facei] < SMALL)
+                {
+                    alphaPhip[facei] = alphap[facei]*phip[facei];
+                }
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H

@@ -319,6 +319,9 @@ public:
             return alphaRhoPhi_;
         }
 
+        //- Ensure that the flux at inflow BCs is preserved
+        void correctInflowFlux(surfaceScalarField& alphaPhi) const;
+
         //- Correct the phase properties
         //  other than the thermodynamics and turbulence
         //  which have special treatment

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -444,28 +444,7 @@ void Foam::twoPhaseSystem::solve()
             )
         );
 
-        surfaceScalarField::Boundary& alphaPhic1Bf =
-            alphaPhic1.boundaryFieldRef();
-
-        // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhic1Bf, patchi)
-        {
-            fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi];
-
-            if (!alphaPhic1p.coupled())
-            {
-                const scalarField& phi1p = phi1.boundaryField()[patchi];
-                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                forAll(alphaPhic1p, facei)
-                {
-                    if (phi1p[facei] < 0)
-                    {
-                        alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
-                    }
-                }
-            }
-        }
+        phase1_.correctInflowFlux(alphaPhic1);
 
         if (nAlphaSubCycles > 1)
         {
@@ -537,6 +516,7 @@ void Foam::twoPhaseSystem::solve()
 
         phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
         alpha2 = scalar(1) - alpha1;
+        phase2_.correctInflowFlux(phase2_.alphaPhi());
         phase2_.alphaRhoPhi() =
             fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();