reactingMultiphaseEulerFoam: New Euler-Euler multiphase solver

Supporting any number of phases with heat and mass transfer, phase-change and reactions

applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H

@@ -18,7 +18,7 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
         new fvVectorMatrix
         (
             fvm::ddt(alpha, U)
-          + fvm::div(phase.phiAlpha(), U)
+          + fvm::div(phase.alphaPhi(), U)
 
           + (alpha/phase.rho())*fluid.Cvm(phase)*
             (

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

@@ -61,7 +61,7 @@ void Foam::multiphaseSystem::calcAlphas()
 
 void Foam::multiphaseSystem::solveAlphas()
 {
-    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
+    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
     int phasei = 0;
 
     forAllIter(PtrDictionary<phaseModel>, phases_, iter)
@@ -69,10 +69,10 @@ void Foam::multiphaseSystem::solveAlphas()
         phaseModel& phase1 = iter();
         volScalarField& alpha1 = phase1;
 
-        phase1.phiAlpha() =
+        phase1.alphaPhi() =
             dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
 
-        phiAlphaCorrs.set
+        alphaPhiCorrs.set
         (
             phasei,
             new surfaceScalarField
@@ -87,7 +87,7 @@ void Foam::multiphaseSystem::solveAlphas()
             )
         );
 
-        surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
+        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
 
         forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
         {
@@ -118,7 +118,7 @@ void Foam::multiphaseSystem::solveAlphas()
                 "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
             );
 
-            phiAlphaCorr += fvc::flux
+            alphaPhiCorr += fvc::flux
             (
                 -fvc::flux(-phir, phase2, phirScheme),
                 phase1,
@@ -127,21 +127,21 @@ void Foam::multiphaseSystem::solveAlphas()
         }
 
         // Ensure that the flux at inflow BCs is preserved
-        forAll(phiAlphaCorr.boundaryField(), patchi)
+        forAll(alphaPhiCorr.boundaryField(), patchi)
         {
-            fvsPatchScalarField& phiAlphaCorrp =
+            fvsPatchScalarField& alphaPhiCorrp =
-                phiAlphaCorr.boundaryField()[patchi];
+                alphaPhiCorr.boundaryField()[patchi];
 
-            if (!phiAlphaCorrp.coupled())
+            if (!alphaPhiCorrp.coupled())
             {
                 const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
                 const scalarField& alpha1p = alpha1.boundaryField()[patchi];
 
-                forAll(phiAlphaCorrp, facei)
+                forAll(alphaPhiCorrp, facei)
                 {
                     if (phi1p[facei] < 0)
                     {
-                        phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei];
+                        alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
                     }
                 }
             }
@@ -153,7 +153,7 @@ void Foam::multiphaseSystem::solveAlphas()
             geometricOneField(),
             phase1,
             phi_,
-            phiAlphaCorr,
+            alphaPhiCorr,
             zeroField(),
             zeroField(),
             1,
@@ -164,7 +164,7 @@ void Foam::multiphaseSystem::solveAlphas()
         phasei++;
     }
 
-    MULES::limitSum(phiAlphaCorrs);
+    MULES::limitSum(alphaPhiCorrs);
 
     volScalarField sumAlpha
     (
@@ -184,19 +184,19 @@ void Foam::multiphaseSystem::solveAlphas()
     {
         phaseModel& phase1 = iter();
 
-        surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
-        phiAlpha += upwind<scalar>(mesh_, phi_).flux(phase1);
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);
 
         MULES::explicitSolve
         (
             geometricOneField(),
             phase1,
-            phiAlpha,
+            alphaPhi,
             zeroField(),
             zeroField()
         );
 
-        phase1.phiAlpha() += phiAlpha;
+        phase1.alphaPhi() += alphaPhi;
 
         Info<< phase1.name() << " volume fraction, min, max = "
             << phase1.weightedAverage(mesh_.V()).value()

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

@@ -100,11 +100,11 @@ Foam::phaseModel::phaseModel
         mesh,
         dimensionedVector("0", dimVelocity/dimTime, vector::zero)
     ),
-    phiAlpha_
+    alphaPhi_
     (
         IOobject
         (
-            IOobject::groupName("phiAlpha", phaseName),
+            IOobject::groupName("alphaPhi", phaseName),
             mesh.time().timeName(),
             mesh
         ),

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.H

@@ -80,7 +80,7 @@ class phaseModel
         volVectorField DDtU_;
 
         //- Volumetric flux of the phase
-        surfaceScalarField phiAlpha_;
+        surfaceScalarField alphaPhi_;
 
         //- Volumetric flux for the phase
         autoPtr<surfaceScalarField> phiPtr_;
@@ -198,14 +198,14 @@ public:
             return phiPtr_();
         }
 
-        const surfaceScalarField& phiAlpha() const
+        const surfaceScalarField& alphaPhi() const
         {
-            return phiAlpha_;
+            return alphaPhi_;
         }
 
-        surfaceScalarField& phiAlpha()
+        surfaceScalarField& alphaPhi()
         {
-            return phiAlpha_;
+            return alphaPhi_;
         }
 
         //- Correct the phase properties

applications/solvers/multiphase/reactingEulerFoam/Allwclean

@@ -6,5 +6,6 @@ wclean libso phaseSystems
 wclean libso interfacialModels
 wclean libso interfacialCompositionModels
 reactingTwoPhaseEulerFoam/Allwclean
+reactingMultiphaseEulerFoam/Allwclean
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/reactingEulerFoam/Allwmake

@@ -8,5 +8,6 @@ wmake libso phaseSystems
 wmake libso interfacialModels
 wmake libso interfacialCompositionModels
 reactingTwoPhaseEulerFoam/Allwmake
+reactingMultiphaseEulerFoam/Allwmake
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/options

@@ -21,7 +21,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
-    -lreactingTwoPhaseSystem \
     -lfluidThermophysicalModels \
     -lreactionThermophysicalModels \
     -lspecie

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/Make/options

@@ -10,7 +10,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
-    -lreactingTwoPhaseSystem \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie

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

@@ -132,7 +132,10 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
     volScalarField muAlphaI
     (
         alpha1*rho1*nu1*alpha2*rho2*nu2
-       /(alpha1*rho1*nu1 + alpha2*rho2*nu2)
+       /(
+           max(alpha1, pair_.phase1().residualAlpha())*rho1*nu1
+         + max(alpha2, pair_.phase2().residualAlpha())*rho2*nu2
+        )
     );
 
     volScalarField ReI

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

@@ -164,6 +164,60 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
 }
 
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::dmdt
+(
+    const Foam::phaseModel& phase
+) const
+{
+    tmp<volScalarField> tdmdt
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("dmdt", phase.name()),
+                this->mesh_.time().timeName(),
+                this->mesh_
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimDensity/dimTime, 0)
+        )
+    );
+
+    forAllConstIter
+    (
+        phaseSystem::phasePairTable,
+        this->phasePairs_,
+        phasePairIter
+    )
+    {
+        const phasePair& pair(phasePairIter());
+
+        if (pair.ordered())
+        {
+            continue;
+        }
+
+        const phaseModel* phase1 = &pair.phase1();
+        const phaseModel* phase2 = &pair.phase2();
+
+        forAllConstIter(phasePair, pair, iter)
+        {
+            if (phase1 == &phase)
+            {
+                tdmdt() += this->dmdt(pair);
+            }
+
+            Swap(phase1, phase2);
+        }
+    }
+
+    return tdmdt;
+}
+
+
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
 Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.H

@@ -121,10 +121,10 @@ public:
     // Member Functions
 
         //- Return the interfacial mass flow rate
-        virtual tmp<volScalarField> dmdt
+        virtual tmp<volScalarField> dmdt(const phasePairKey& key) const;
-        (
+
-            const phasePairKey& key
+        //- Return the total interfacial mass transfer rate for phase
-        ) const;
+        virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
 
         //- Return the momentum transfer matrices
         virtual autoPtr<phaseSystem::momentumTransferTable>

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

@@ -85,6 +85,30 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
 }
 
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::HeatTransferPhaseSystem<BasePhaseSystem>::dmdt
+(
+    const Foam::phaseModel& phase
+) const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("dmdt", phase.name()),
+                this->mesh().time().timeName(),
+                this->mesh().time()
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimDensity/dimTime, 0)
+        )
+    );
+}
+
+
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::heatTransferTable>
 Foam::HeatTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.H

@@ -95,6 +95,9 @@ public:
             const phasePairKey& key
         ) const;
 
+        //- Return the total interfacial mass transfer rate for phase
+        virtual tmp<volScalarField> dmdt(const phaseModel& phase) const;
+
         //- Return the heat transfer matrices
         virtual autoPtr<phaseSystem::heatTransferTable>
             heatTransfer() const;

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

@@ -182,7 +182,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kd
         )
     );
 
-    // Add the implicit part of the drag force
     forAllConstIter
     (
         phaseSystem::KdTable,
@@ -442,7 +441,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
     }
 
     // Add the implicit part of the drag force
-    /* ***HGW Currently this is handled in the pEqn
     forAllConstIter
     (
         phaseSystem::KdTable,
@@ -466,7 +464,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
             Swap(phase, otherPhase);
         }
     }
-    */
 
     // Update the virtual mass coefficients
     forAllConstIter
@@ -499,20 +496,57 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
             const volVectorField& U = phase->U();
             const surfaceScalarField& phi = phase->phi();
 
-            *eqns[phase->name()] +=
+            *eqns[phase->name()] -=
-              - Vm
+                Vm
                *(
                     fvm::ddt(U)
                   + fvm::div(phi, U)
                   - fvm::Sp(fvc::div(phi), U)
                   - otherPhase->DUDt()
                 )
-              - this->MRF_.DDt(Vm, U - otherPhase->U());
+              + this->MRF_.DDt(Vm, U - otherPhase->U());
 
             Swap(phase, otherPhase);
         }
     }
 
+    return eqnsPtr;
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::PtrList<Foam::volVectorField> >
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
+{
+    autoPtr<PtrList<volVectorField> > tFs
+    (
+        new PtrList<volVectorField>(this->phases().size())
+    );
+
+    PtrList<volVectorField>& Fs = tFs();
+
+    forAll(Fs, phasei)
+    {
+        Fs.set
+        (
+            phasei,
+            new volVectorField
+            (
+                IOobject
+                (
+                    liftModel::typeName + ":F",
+                    this->mesh_.time().timeName(),
+                    this->mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                this->mesh_,
+                dimensionedVector("zero", liftModel::dimF, vector::zero)
+            )
+        );
+    }
+
     // Add the lift force
     forAllConstIter
     (
@@ -525,8 +559,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 
         const phasePair& pair(this->phasePairs_[liftModelIter.key()]);
 
-        *eqns[pair.phase1().name()] += F;
+        Fs[pair.phase1().index()] += F;
-        *eqns[pair.phase2().name()] -= F;
+        Fs[pair.phase2().index()] -= F;
     }
 
     // Add the wall lubrication force
@@ -542,28 +576,28 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
         const phasePair&
             pair(this->phasePairs_[wallLubricationModelIter.key()]);
 
-        *eqns[pair.phase1().name()] += F;
+        Fs[pair.phase1().index()] += F;
-        *eqns[pair.phase2().name()] -= F;
+        Fs[pair.phase2().index()] -= F;
     }
 
     // Add the turbulent dispersion force
-    forAllConstIter
+    // forAllConstIter
-    (
+    // (
-        turbulentDispersionModelTable,
+    //     turbulentDispersionModelTable,
-        turbulentDispersionModels_,
+    //     turbulentDispersionModels_,
-        turbulentDispersionModelIter
+    //     turbulentDispersionModelIter
-    )
+    // )
-    {
+    // {
-        const volVectorField F(turbulentDispersionModelIter()->F<vector>());
+    //     const volVectorField F(turbulentDispersionModelIter()->F<vector>());
 
-        const phasePair&
+    //     const phasePair&
-            pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
+    //         pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
 
-        *eqns[pair.phase1().name()] += F;
+    //     *eqns[pair.phase1().name()] += F;
-        *eqns[pair.phase2().name()] -= F;
+    //     *eqns[pair.phase2().name()] -= F;
-    }
+    // }
 
-    return eqnsPtr;
+    return tFs;
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H

@@ -172,6 +172,9 @@ public:
         //- Return the combined force (lift + wall-lubrication)
         virtual tmp<volVectorField> F(const phasePairKey& key) const;
 
+        //- Return the combined force (lift + wall-lubrication)
+        virtual autoPtr<PtrList<volVectorField> > Fs() const;
+
         //- Return the combined face-force (lift + wall-lubrication)
         virtual tmp<surfaceScalarField> Ff(const phasePairKey& key) const;
 

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

@@ -32,10 +32,11 @@ template<class BasePhaseModel>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
-    BasePhaseModel(fluid, phaseName),
+    BasePhaseModel(fluid, phaseName, index),
     divU_
     (
         IOobject
@@ -134,7 +135,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 template<class BasePhaseModel>
 bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const
 {
-    return true;
+    return !this->thermo().incompressible();
 }
 
 

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

@@ -65,7 +65,12 @@ public:
 
     // Constructors
 
-        AnisothermalPhaseModel(const phaseSystem& fluid, const word& phaseName);
+        AnisothermalPhaseModel
+        (
+            const phaseSystem& fluid,
+            const word& phaseName,
+            const label index
+        );
 
 
     //- Destructor

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

@@ -32,10 +32,11 @@ template<class BasePhaseModel>
 Foam::InertPhaseModel<BasePhaseModel>::InertPhaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
-    BasePhaseModel(fluid, phaseName)
+    BasePhaseModel(fluid, phaseName, index)
 {}
 
 

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

@@ -56,7 +56,12 @@ public:
 
     // Constructors
 
-        InertPhaseModel(const phaseSystem& fluid, const word& phaseName);
+        InertPhaseModel
+        (
+            const phaseSystem& fluid,
+            const word& phaseName,
+            const label index
+        );
 
 
     //- Destructor

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

@@ -32,10 +32,11 @@ template<class BasePhaseModel>
 Foam::IsothermalPhaseModel<BasePhaseModel>::IsothermalPhaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
-    BasePhaseModel(fluid, phaseName)
+    BasePhaseModel(fluid, phaseName, index)
 {}
 
 

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

@@ -57,7 +57,12 @@ public:
 
     // Constructors
 
-        IsothermalPhaseModel(const phaseSystem& fluid, const word& phaseName);
+        IsothermalPhaseModel
+        (
+            const phaseSystem& fluid,
+            const word& phaseName,
+            const label index
+        );
 
 
     //- Destructor

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

@@ -123,10 +123,11 @@ template<class BasePhaseModel>
 Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
-    BasePhaseModel(fluid, phaseName),
+    BasePhaseModel(fluid, phaseName, index),
     U_
     (
         IOobject
@@ -197,6 +198,7 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
         dimensionedScalar("0", dimDensity/dimTime, 0)
     )
 {
+    phi_.writeOpt() = IOobject::AUTO_WRITE;
     correctKinematics();
 }
 

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

@@ -106,7 +106,8 @@ public:
         MovingPhaseModel
         (
             const phaseSystem& fluid,
-            const word& phaseName
+            const word& phaseName,
+            const label index
         );
 
 

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

@@ -41,10 +41,11 @@ template<class BasePhaseModel>
 Foam::MultiComponentPhaseModel<BasePhaseModel>::MultiComponentPhaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
-    BasePhaseModel(fluid, phaseName),
+    BasePhaseModel(fluid, phaseName, index),
     Sc_
     (
         "Sc",

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

@@ -73,7 +73,8 @@ public:
         MultiComponentPhaseModel
         (
             const phaseSystem& fluid,
-            const word& phaseName
+            const word& phaseName,
+            const label index
         );
 
 

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

@@ -32,10 +32,11 @@ template<class BasePhaseModel>
 Foam::PurePhaseModel<BasePhaseModel>::PurePhaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
-    BasePhaseModel(fluid, phaseName)
+    BasePhaseModel(fluid, phaseName, index)
 {}
 
 

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

@@ -64,7 +64,12 @@ public:
 
     // Constructors
 
-        PurePhaseModel(const phaseSystem& fluid, const word& phaseName);
+        PurePhaseModel
+        (
+            const phaseSystem& fluid,
+            const word& phaseName,
+            const label index
+        );
 
 
     //- Destructor

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

@@ -33,10 +33,11 @@ template<class BasePhaseModel, class ReactionType>
 Foam::ReactingPhaseModel<BasePhaseModel, ReactionType>::ReactingPhaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
-    BasePhaseModel(fluid, phaseName, false),
+    BasePhaseModel(fluid, phaseName, index, false),
     reaction_
     (
         ReactionType::New(fluid.mesh(), this->name())

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

@@ -67,7 +67,8 @@ public:
         ReactingPhaseModel
         (
             const phaseSystem& fluid,
-            const word& phaseName
+            const word& phaseName,
+            const label index
         );
 
 

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

@@ -41,10 +41,11 @@ Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::ThermoPhaseModel
 (
     const phaseSystem& fluid,
     const word& phaseName,
+    const label index,
     const bool createThermo
 )
 :
-    BasePhaseModel(fluid, phaseName)
+    BasePhaseModel(fluid, phaseName, index)
 {
     if (createThermo)
     {

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

@@ -76,6 +76,7 @@ public:
         (
             const phaseSystem& fluid,
             const word& phaseName,
+            const label index,
             const bool createThermo = true
         );
 

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

@@ -31,7 +31,8 @@ License
 Foam::autoPtr<Foam::phaseModel> Foam::phaseModel::New
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 {
     word phaseModelType(fluid.subDict(phaseName).lookup("type"));
@@ -52,7 +53,7 @@ Foam::autoPtr<Foam::phaseModel> Foam::phaseModel::New
             << exit(FatalError);
     }
 
-    return cstrIter()(fluid, phaseName);
+    return cstrIter()(fluid, phaseName, index);
 }
 
 // ************************************************************************* //

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

@@ -41,7 +41,8 @@ namespace Foam
 Foam::phaseModel::phaseModel
 (
     const phaseSystem& fluid,
-    const word& phaseName
+    const word& phaseName,
+    const label index
 )
 :
     volScalarField
@@ -60,6 +61,7 @@ Foam::phaseModel::phaseModel
 
     fluid_(fluid),
     name_(phaseName),
+    index_(index),
     residualAlpha_
     (
         "residualAlpha",
@@ -99,6 +101,12 @@ const Foam::word& Foam::phaseModel::keyword() const
 }
 
 
+Foam::label Foam::phaseModel::index() const
+{
+    return index_;
+}
+
+
 const Foam::phaseSystem& Foam::phaseModel::fluid() const
 {
     return fluid_;

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

@@ -67,6 +67,9 @@ class phaseModel
         //- Name of phase
         word name_;
 
+        //- Index of phase
+        label index_;
+
         //- Return the residual phase-fraction for given phase
         //  Used to stabilize the phase momentum as the phase-fraction -> 0
         dimensionedScalar residualAlpha_;
@@ -93,9 +96,10 @@ public:
             phaseSystem,
             (
                 const phaseSystem& fluid,
-                const word& phaseName
+                const word& phaseName,
+                const label index
             ),
-            (fluid, phaseName)
+            (fluid, phaseName, index)
         );
 
 
@@ -104,7 +108,8 @@ public:
         phaseModel
         (
             const phaseSystem& fluid,
-            const word& phaseName
+            const word& phaseName,
+            const label index
         );
 
         //- Return clone
@@ -116,7 +121,8 @@ public:
         static autoPtr<phaseModel> New
         (
             const phaseSystem& fluid,
-            const word& phaseName
+            const word& phaseName,
+            const label index
         );
 
         //- Return a pointer to a new phase created on freestore
@@ -124,6 +130,7 @@ public:
         class iNew
         {
             const phaseSystem& fluid_;
+            mutable label indexCounter_;
 
         public:
 
@@ -132,14 +139,16 @@ public:
                 const phaseSystem& fluid
             )
             :
-                fluid_(fluid)
+                fluid_(fluid),
+                indexCounter_(-1)
             {}
 
             autoPtr<phaseModel> operator()(Istream& is) const
             {
+                indexCounter_++;
                 return autoPtr<phaseModel>
                 (
-                    phaseModel::New(fluid_, word(is))
+                    phaseModel::New(fluid_, word(is), indexCounter_)
                 );
             }
         };
@@ -154,8 +163,12 @@ public:
         //- Return the name of this phase
         const word& name() const;
 
+        //- Return the name of the phase for use as the keyword in PtrDictionary
         const word& keyword() const;
 
+        //- Return the index of the phase
+        label index() const;
+
         //- Return the system to which this phase belongs
         const phaseSystem& fluid() const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C

@@ -162,6 +162,8 @@ Foam::phaseSystem::phaseSystem
     MRF_(mesh_),
     fvOptions_(mesh_)
 {
+    phi_.writeOpt() = IOobject::AUTO_WRITE;
+
     // Blending methods
     forAllConstIter(dictionary, subDict("blending"), iter)
     {

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Allwclean

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

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Allwmake

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wmake libso multiphaseSystem
+wmake libso multiphaseCompressibleTurbulenceModels
+wmake
+
+# ----------------------------------------------------------------- end-of-file

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

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar CoNum = 0.0;
+scalar meanCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        fvc::surfaceSum(mag(phi))().internalField()
+    );
+
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        sumPhi = max
+        (
+            sumPhi,
+            fvc::surfaceSum(mag(iter().phi()))().internalField()
+        );
+    }
+
+    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Courant Number mean: " << meanCoNum
+    << " max: " << CoNum << endl;
+
+// ************************************************************************* //

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

@@ -0,0 +1,47 @@
+{
+    fluid.correctEnergyTransport();
+
+    autoPtr<phaseSystem::heatTransferTable>
+        heatTransferPtr(fluid.heatTransfer());
+
+    phaseSystem::heatTransferTable& heatTransfer = heatTransferPtr();
+
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        const volScalarField& alpha = phase;
+        const volScalarField& rho = phase.rho();
+        const volVectorField& U = phase.U();
+
+        tmp<fvScalarMatrix> EEqn(phase.heEqn());
+
+        if (EEqn.valid())
+        {
+            EEqn =
+            (
+                EEqn
+             ==
+               *heatTransfer[phase.name()]
+              + alpha*rho*(U&g)
+              + fvOptions(alpha, rho, phase.thermo().he())
+            );
+
+            EEqn->relax();
+            fvOptions.constrain(EEqn());
+            EEqn->solve();
+        }
+    }
+}
+
+fluid.correctThermo();
+
+forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+{
+    phaseModel& phase = iter();
+
+    Info<< phase.name() << " min/max T "
+        << min(phase.thermo().T()).value()
+        << " - "
+        << max(phase.thermo().T()).value()
+        << endl;
+}

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/files

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

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options

@@ -0,0 +1,25 @@
+EXE_INC = -ggdb3 \
+    -ImultiphaseSystem/lnInclude \
+    -I../phaseSystems/lnInclude \
+    -I../interfacialModels/lnInclude \
+    -I../interfacialCompositionModels/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -lreactingPhaseSystem \
+    -lreactingMultiphaseSystem \
+    -lreactingEulerianInterfacialModels \
+    -lreactingEulerianInterfacialCompositionModels \
+    -lmultiphaseReactingTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling

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

@@ -0,0 +1,36 @@
+{
+    autoPtr<phaseSystem::massTransferTable>
+        massTransferPtr(fluid.massTransfer());
+
+    phaseSystem::massTransferTable&
+        massTransfer(massTransferPtr());
+
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        PtrList<volScalarField>& Y = phase.Y();
+        const volScalarField& alpha = phase;
+        const volScalarField& rho = phase.rho();
+
+        forAll(Y, i)
+        {
+            tmp<fvScalarMatrix> YiEqn(phase.YiEqn(Y[i]));
+
+            if (YiEqn.valid())
+            {
+                YiEqn =
+                (
+                    YiEqn
+                  ==
+                    *massTransfer[Y[i].name()]
+                   + fvOptions(alpha, rho, Y[i])
+                );
+
+                YiEqn->relax();
+                YiEqn->solve(mesh.solver("Yi"));
+            }
+        }
+    }
+
+    fluid.massTransfer(); // updates interfacial mass flow rates
+}

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

@@ -0,0 +1,52 @@
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+
+Info<< "Creating phaseSystem\n" << endl;
+
+autoPtr<multiphaseSystem> fluidPtr
+(
+    multiphaseSystem::New(mesh)
+);
+multiphaseSystem& fluid = fluidPtr();
+
+surfaceScalarField& phi = fluid.phi();
+
+dimensionedScalar pMin
+(
+    "pMin",
+    dimPressure,
+    fluid
+);
+
+#include "gh.H"
+
+volScalarField& p = fluid.phases().first().thermo().p();
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+mesh.setFluxRequired(p_rgh.name());
+
+const IOMRFZoneList& MRF = fluid.MRF();
+fv::IOoptionList& fvOptions = fluid.fvOptions();

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/log

@@ -0,0 +1,6 @@
+/home/dm2/henry/OpenFOAM/OpenFOAM-dev/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam
+Making dependency list for source file reactingMultiphaseEulerFoam.C
+clang++ -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-c++11-extensions -O3  -DNoRepository -ftemplate-depth-100 -ImultiphaseSystem/lnInclude -I../phaseSystems/lnInclude -I../interfacialModels/lnInclude -I../interfacialCompositionModels/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/basic/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/transportModels/compressible/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/TurbulenceModels/turbulenceModels/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/TurbulenceModels/compressible/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/TurbulenceModels/phaseCompressible/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/finiteVolume/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/fvOptions/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/meshTools/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/sampling/lnInclude -IlnInclude -I. -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/OpenFOAM/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/OSspecific/POSIX/lnInclude   -fPIC -c reactingMultiphaseEulerFoam.C -o /home/dm2/henry/OpenFOAM/OpenFOAM-dev/platforms/linux64ClangDPInt32Opt/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.o
+clang++ -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-c++11-extensions -O3  -DNoRepository -ftemplate-depth-100 -ImultiphaseSystem/lnInclude -I../phaseSystems/lnInclude -I../interfacialModels/lnInclude -I../interfacialCompositionModels/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/basic/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/transportModels/compressible/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/TurbulenceModels/turbulenceModels/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/TurbulenceModels/compressible/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/TurbulenceModels/phaseCompressible/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/finiteVolume/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/fvOptions/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/meshTools/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/sampling/lnInclude -IlnInclude -I. -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/OpenFOAM/lnInclude -I/home/dm2/henry/OpenFOAM/OpenFOAM-dev/src/OSspecific/POSIX/lnInclude   -fPIC -Xlinker --add-needed /home/dm2/henry/OpenFOAM/OpenFOAM-dev/platforms/linux64ClangDPInt32Opt/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.o -L/home/dm2/henry/OpenFOAM/OpenFOAM-dev/platforms/linux64ClangDPInt32Opt/lib \
+    -lreactingPhaseSystem -lreactingMultiphaseSystem -lreactingEulerianInterfacialModels -lreactingEulerianInterfacialCompositionModels -lmultiphaseReactingTurbulenceModels -lfiniteVolume -lfvOptions -lmeshTools -lsampling -lOpenFOAM -ldl  \
+     -lm -o /home/dm2/henry/OpenFOAM/OpenFOAM-dev/platforms/linux64ClangDPInt32Opt/bin/reactingMultiphaseEulerFoam

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/Make/files

@@ -0,0 +1,3 @@
+multiphaseCompressibleTurbulenceModels.C
+
+LIB = $(FOAM_LIBBIN)/libmultiphaseReactingTurbulenceModels

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/Make/options

@@ -0,0 +1,13 @@
+EXE_INC = \
+    -I../multiphaseSystem/lnInclude \
+    -I../../phaseSystems/lnInclude \
+    -I../../interfacialModels/lnInclude\
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/transportModel \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude

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

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014-2015 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 "PhaseCompressibleTurbulenceModel.H"
+#include "phaseModel.H"
+#include "multiphaseSystem.H"
+#include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
+
+#include "ThermalDiffusivity.H"
+#include "EddyDiffusivity.H"
+
+#include "laminar.H"
+#include "RASModel.H"
+#include "LESModel.H"
+
+makeBaseTurbulenceModel
+(
+    volScalarField,
+    volScalarField,
+    compressibleTurbulenceModel,
+    PhaseCompressibleTurbulenceModel,
+    ThermalDiffusivity,
+    phaseModel
+);
+
+#define makeRASModel(Type)                                                     \
+    makeTemplatedTurbulenceModel                                               \
+    (phaseModelPhaseCompressibleTurbulenceModel, RAS, Type)
+
+#define makeLESModel(Type)                                                     \
+    makeTemplatedTurbulenceModel                                               \
+    (phaseModelPhaseCompressibleTurbulenceModel, LES, Type)
+
+#include "kEpsilon.H"
+makeRASModel(kEpsilon);
+
+#include "kOmegaSST.H"
+makeRASModel(kOmegaSST);
+
+#include "Smagorinsky.H"
+makeLESModel(Smagorinsky);
+
+#include "kEqn.H"
+makeLESModel(kEqn);
+
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/Make/files

@@ -0,0 +1,6 @@
+alphaContactAngle/alphaContactAngleFvPatchScalarField.C
+multiphaseSystem.C
+newMultiphaseSystem.C
+multiphaseSystems.C
+
+LIB = $(FOAM_LIBBIN)/libreactingMultiphaseSystem

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/Make/options

@@ -0,0 +1,21 @@
+EXE_INC = -ggdb3 \
+    -I../../interfacialModels/lnInclude \
+    -I../../interfacialCompositionModels/lnInclude \
+    -I../../phaseSystems/lnInclude \
+    -IalphaContactAngle \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lcombustionModels \
+    -lreactingPhaseSystem

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/alphaContactAngle/alphaContactAngleFvPatchScalarField.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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 "alphaContactAngleFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
+(
+    Istream& is
+)
+:
+    theta0_(readScalar(is)),
+    uTheta_(readScalar(is)),
+    thetaA_(readScalar(is)),
+    thetaR_(readScalar(is))
+{}
+
+
+Istream& operator>>
+(
+    Istream& is,
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    return is;
+}
+
+
+Ostream& operator<<
+(
+    Ostream& os,
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    os  << tp.theta0_ << token::SPACE
+        << tp.uTheta_ << token::SPACE
+        << tp.thetaA_ << token::SPACE
+        << tp.thetaR_;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(p, iF)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    thetaProps_(dict.lookup("thetaProperties"))
+{
+    evaluate();
+}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, iF),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("thetaProperties")
+        << thetaProps_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphaContactAngleFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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::alphaContactAngleFvPatchScalarField
+
+Description
+    Contact-angle boundary condition for multi-phase interface-capturing
+    simulations.  Used in conjuction with multiphaseSystem.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class alphaContactAngleFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphaContactAngleFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+public:
+
+    class interfaceThetaProps
+    {
+        //- Equilibrium contact angle
+        scalar theta0_;
+
+        //- Dynamic contact angle velocity scale
+        scalar uTheta_;
+
+        //- Limiting advancing contact angle
+        scalar thetaA_;
+
+        //- Limiting receeding contact angle
+        scalar thetaR_;
+
+
+    public:
+
+        // Constructors
+            interfaceThetaProps()
+            {}
+
+            interfaceThetaProps(Istream&);
+
+
+        // Member functions
+
+            //- Return the equilibrium contact angle theta0
+            scalar theta0(bool matched=true) const
+            {
+                if (matched) return theta0_;
+                else return 180.0 - theta0_;
+            }
+
+            //- Return the dynamic contact angle velocity scale
+            scalar uTheta() const
+            {
+                return uTheta_;
+            }
+
+            //- Return the limiting advancing contact angle
+            scalar thetaA(bool matched=true) const
+            {
+                if (matched) return thetaA_;
+                else return 180.0 - thetaA_;
+            }
+
+            //- Return the limiting receeding contact angle
+            scalar thetaR(bool matched=true) const
+            {
+                if (matched) return thetaR_;
+                else return 180.0 - thetaR_;
+            }
+
+
+        // IO functions
+
+            friend Istream& operator>>(Istream&, interfaceThetaProps&);
+            friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
+    };
+
+    typedef HashTable
+    <
+        interfaceThetaProps,
+        phasePairKey,
+        phasePairKey::hash
+    > thetaPropsTable;
+
+
+private:
+
+    // Private data
+
+        thetaPropsTable thetaProps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphaContactAngle");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
+        //  onto a new patch
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return the contact angle properties
+        const thetaPropsTable& thetaProps() const
+        {
+            return thetaProps_;
+        }
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,713 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2015 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 "multiphaseSystem.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+
+#include "MULES.H"
+#include "subCycle.H"
+
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvcSnGrad.H"
+#include "fvcFlux.H"
+#include "fvcMeshPhi.H"
+#include "fvcSup.H"
+
+#include "fvmDdt.H"
+#include "fvmLaplacian.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(multiphaseSystem, 0);
+    defineRunTimeSelectionTable(multiphaseSystem, dictionary);
+}
+
+const Foam::scalar Foam::multiphaseSystem::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::multiphaseSystem::calcAlphas()
+{
+    scalar level = 0.0;
+    alphas_ == 0.0;
+
+    forAllIter(PtrDictionary<phaseModel>, phases(), iter)
+    {
+        alphas_ += level*iter();
+        level += 1.0;
+    }
+
+    alphas_.correctBoundaryConditions();
+}
+
+
+void Foam::multiphaseSystem::solveAlphas()
+{
+    PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
+    int phasei = 0;
+
+    forAllIter(PtrDictionary<phaseModel>, phases(), iter)
+    {
+        phaseModel& phase = iter();
+        volScalarField& alpha1 = phase;
+
+        phase.alphaPhi() =
+            dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
+
+        alphaPhiCorrs.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                "phi" + alpha1.name() + "Corr",
+                fvc::flux
+                (
+                    phi_,
+                    phase,
+                    "div(phi," + alpha1.name() + ')'
+                )
+            )
+        );
+
+        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
+
+        forAllIter(PtrDictionary<phaseModel>, phases(), iter2)
+        {
+            phaseModel& phase2 = iter2();
+            volScalarField& alpha2 = phase2;
+
+            if (&phase2 == &phase) continue;
+
+            surfaceScalarField phir(phase.phi() - phase2.phi());
+
+            cAlphaTable::const_iterator cAlpha
+            (
+                cAlphas_.find(phasePairKey(phase.name(), phase2.name()))
+            );
+
+            if (cAlpha != cAlphas_.end())
+            {
+                surfaceScalarField phic
+                (
+                    (mag(phi_) + mag(phir))/mesh_.magSf()
+                );
+
+                phir += min(cAlpha()*phic, max(phic))*nHatf(phase, phase2);
+            }
+
+            word phirScheme
+            (
+                "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
+            );
+
+            alphaPhiCorr += fvc::flux
+            (
+                -fvc::flux(-phir, phase2, phirScheme),
+                phase,
+                phirScheme
+            );
+        }
+
+        // Ensure that the flux at inflow BCs is preserved
+        forAll(alphaPhiCorr.boundaryField(), patchi)
+        {
+            fvsPatchScalarField& alphaPhiCorrp =
+                alphaPhiCorr.boundaryField()[patchi];
+
+            if (!alphaPhiCorrp.coupled())
+            {
+                const scalarField& phi1p = phase.phi().boundaryField()[patchi];
+                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
+
+                forAll(alphaPhiCorrp, facei)
+                {
+                    if (phi1p[facei] < 0)
+                    {
+                        alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
+                    }
+                }
+            }
+        }
+
+        MULES::limit
+        (
+            1.0/mesh_.time().deltaT().value(),
+            geometricOneField(),
+            phase,
+            phi_,
+            alphaPhiCorr,
+            zeroField(),
+            zeroField(),
+            phase.alphaMax(),
+            0,
+            true
+        );
+
+        phasei++;
+    }
+
+    MULES::limitSum(alphaPhiCorrs);
+
+    volScalarField sumAlpha
+    (
+        IOobject
+        (
+            "sumAlpha",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedScalar("sumAlpha", dimless, 0)
+    );
+
+
+    volScalarField divU(fvc::div(fvc::absolute(phi_, phases().first().U())));
+
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, phases(), iter)
+    {
+        phaseModel& phase = iter();
+        volScalarField& alpha = phase;
+
+        surfaceScalarField& alphaPhic = alphaPhiCorrs[phasei];
+        alphaPhic += upwind<scalar>(mesh_, phi_).flux(phase);
+
+        volScalarField::DimensionedInternalField Sp
+        (
+            IOobject
+            (
+                "Sp",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("Sp", phase.divU().dimensions(), 0.0)
+        );
+
+        volScalarField::DimensionedInternalField Su
+        (
+            IOobject
+            (
+                "Su",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            // Divergence term is handled explicitly to be
+            // consistent with the explicit transport solution
+            divU*min(alpha, scalar(1))
+        );
+
+        if (phase.compressible())
+        {
+            const scalarField& dgdt = phase.divU();
+
+            forAll(dgdt, celli)
+            {
+                if (dgdt[celli] > 0.0)
+                {
+                    Sp[celli] -= dgdt[celli];
+                    Su[celli] += dgdt[celli];
+                }
+                else if (dgdt[celli] < 0.0)
+                {
+                    Sp[celli] +=
+                        dgdt[celli]
+                       *(1.0 - alpha[celli])/max(alpha[celli], 1e-4);
+                }
+            }
+        }
+
+        forAllConstIter(PtrDictionary<phaseModel>, phases(), iter2)
+        {
+            const phaseModel& phase2 = iter2();
+            const volScalarField& alpha2 = phase2;
+
+            if (&phase2 == &phase) continue;
+
+            if (phase2.compressible())
+            {
+                const scalarField& dgdt2 = phase2.divU();
+
+                forAll(dgdt2, celli)
+                {
+                    if (dgdt2[celli] < 0.0)
+                    {
+                        Sp[celli] +=
+                            dgdt2[celli]
+                           *(1.0 - alpha2[celli])/max(alpha2[celli], 1e-4);
+
+                        Su[celli] -=
+                            dgdt2[celli]
+                           *alpha[celli]/max(alpha2[celli], 1e-4);
+                    }
+                    else if (dgdt2[celli] > 0.0)
+                    {
+                        Sp[celli] -= dgdt2[celli];
+                    }
+                }
+            }
+        }
+
+        MULES::explicitSolve
+        (
+            geometricOneField(),
+            alpha,
+            alphaPhic,
+            Sp,
+            Su
+        );
+
+        phase.alphaPhi() += alphaPhic;
+
+        Info<< phase.name() << " volume fraction, min, max = "
+            << phase.weightedAverage(mesh_.V()).value()
+            << ' ' << min(phase).value()
+            << ' ' << max(phase).value()
+            << endl;
+
+        sumAlpha += phase;
+
+        phasei++;
+    }
+
+    Info<< "Phase-sum volume fraction, min, max = "
+        << sumAlpha.weightedAverage(mesh_.V()).value()
+        << ' ' << min(sumAlpha).value()
+        << ' ' << max(sumAlpha).value()
+        << endl;
+}
+
+
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseSystem::nHatfv
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    /*
+    // Cell gradient of alpha
+    volVectorField gradAlpha =
+        alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
+    */
+
+    surfaceVectorField gradAlphaf
+    (
+        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
+      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
+    );
+
+    // Face unit interface normal
+    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::nHatf
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    // Face unit interface normal flux
+    return nHatfv(alpha1, alpha2) & mesh_.Sf();
+}
+
+
+// Correction for the boundary condition on the unit normal nHat on
+// walls to produce the correct contact angle.
+
+// The dynamic contact angle is calculated from the component of the
+// velocity on the direction of the interface, parallel to the wall.
+
+void Foam::multiphaseSystem::correctContactAngle
+(
+    const phaseModel& phase1,
+    const phaseModel& phase2,
+    surfaceVectorField::GeometricBoundaryField& nHatb
+) const
+{
+    const volScalarField::GeometricBoundaryField& gbf
+        = phase1.boundaryField();
+
+    const fvBoundaryMesh& boundary = mesh_.boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& acap =
+                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
+
+            vectorField& nHatPatch = nHatb[patchi];
+
+            vectorField AfHatPatch
+            (
+                mesh_.Sf().boundaryField()[patchi]
+               /mesh_.magSf().boundaryField()[patchi]
+            );
+
+            alphaContactAngleFvPatchScalarField::thetaPropsTable::
+                const_iterator tp =
+                    acap.thetaProps()
+                   .find(phasePairKey(phase1.name(), phase2.name()));
+
+            if (tp == acap.thetaProps().end())
+            {
+                FatalErrorIn
+                (
+                    "multiphaseSystem::correctContactAngle"
+                    "(const phaseModel& phase1, const phaseModel& phase2, "
+                    "fvPatchVectorFieldField& nHatb) const"
+                )   << "Cannot find interface "
+                    << phasePairKey(phase1.name(), phase2.name())
+                    << "\n    in table of theta properties for patch "
+                    << acap.patch().name()
+                    << exit(FatalError);
+            }
+
+            bool matched = (tp.key().first() == phase1.name());
+
+            scalar theta0 = convertToRad*tp().theta0(matched);
+            scalarField theta(boundary[patchi].size(), theta0);
+
+            scalar uTheta = tp().uTheta();
+
+            // Calculate the dynamic contact angle if required
+            if (uTheta > SMALL)
+            {
+                scalar thetaA = convertToRad*tp().thetaA(matched);
+                scalar thetaR = convertToRad*tp().thetaR(matched);
+
+                // Calculated the component of the velocity parallel to the wall
+                vectorField Uwall
+                (
+                    phase1.U()().boundaryField()[patchi].patchInternalField()
+                  - phase1.U()().boundaryField()[patchi]
+                );
+                Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
+
+                // Find the direction of the interface parallel to the wall
+                vectorField nWall
+                (
+                    nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
+                );
+
+                // Normalise nWall
+                nWall /= (mag(nWall) + SMALL);
+
+                // Calculate Uwall resolved normal to the interface parallel to
+                // the interface
+                scalarField uwall(nWall & Uwall);
+
+                theta += (thetaA - thetaR)*tanh(uwall/uTheta);
+            }
+
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12(nHatPatch & AfHatPatch);
+
+            scalarField b1(cos(theta));
+
+            scalarField b2(nHatPatch.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det(1.0 - a12*a12);
+
+            scalarField a((b1 - a12*b2)/det);
+            scalarField b((b2 - a12*b1)/det);
+
+            nHatPatch = a*AfHatPatch + b*nHatPatch;
+
+            nHatPatch /= (mag(nHatPatch) + deltaN_.value());
+        }
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
+(
+    const phaseModel& phase1,
+    const phaseModel& phase2
+) const
+{
+    tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
+
+    correctContactAngle(phase1, phase2, tnHatfv().boundaryField());
+
+    // Simple expression for curvature
+    return -fvc::div(tnHatfv & mesh_.Sf());
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseSystem::multiphaseSystem
+(
+    const fvMesh& mesh
+)
+:
+    phaseSystem(mesh),
+
+    alphas_
+    (
+        IOobject
+        (
+            "alphas",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("alphas", dimless, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+
+    cAlphas_(lookup("interfaceCompression")),
+
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh.V()), 1.0/3.0)
+    )
+{
+    forAllIter(phaseSystem::phaseModelTable, phases(), iter)
+    {
+        volScalarField& alphai = iter();
+        mesh.setFluxRequired(alphai.name());
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::multiphaseSystem::~multiphaseSystem()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
+(
+    const phaseModel& phase1
+) const
+{
+    tmp<surfaceScalarField> tSurfaceTension
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "surfaceTension",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "surfaceTension",
+                dimensionSet(1, -2, -2, 0, 0),
+                0
+            )
+        )
+    );
+
+    forAllConstIter(PtrDictionary<phaseModel>, phases(), iter)
+    {
+        const phaseModel& phase2 = iter();
+
+        if (&phase2 != &phase1)
+        {
+            phasePairKey key12(phase1.name(), phase2.name());
+
+            cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
+
+            if (cAlpha != cAlphas_.end())
+            {
+                tSurfaceTension() +=
+                    fvc::interpolate(sigma(key12)*K(phase1, phase2))
+                   *(
+                        fvc::interpolate(phase2)*fvc::snGrad(phase1)
+                      - fvc::interpolate(phase1)*fvc::snGrad(phase2)
+                    );
+            }
+        }
+    }
+
+    return tSurfaceTension;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseSystem::nearInterface() const
+{
+    tmp<volScalarField> tnearInt
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "nearInterface",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("nearInterface", dimless, 0.0)
+        )
+    );
+
+    forAllConstIter(PtrDictionary<phaseModel>, phases(), iter)
+    {
+        tnearInt() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
+    }
+
+    return tnearInt;
+}
+
+
+void Foam::multiphaseSystem::solve()
+{
+    const fvMesh& mesh = this->mesh();
+    const Time& runTime = mesh.time();
+
+    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+
+    bool LTS = fv::localEulerDdt::enabled(mesh);
+
+    if (nAlphaSubCycles > 1)
+    {
+        tmp<volScalarField> trSubDeltaT;
+
+        if (LTS)
+        {
+            trSubDeltaT =
+                fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+        }
+
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+
+        PtrList<volScalarField> alpha0s(phases().size());
+        PtrList<surfaceScalarField> phiSums(phases().size());
+
+        int phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, phases(), iter)
+        {
+            phaseModel& phase = iter();
+            volScalarField& alpha = phase;
+
+            alpha0s.set
+            (
+                phasei,
+                new volScalarField(alpha.oldTime())
+            );
+
+            phiSums.set
+            (
+                phasei,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "phiSum" + alpha.name(),
+                        runTime.timeName(),
+                        mesh_
+                    ),
+                    mesh_,
+                    dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
+                )
+            );
+
+            phasei++;
+        }
+
+        for
+        (
+            subCycleTime alphaSubCycle
+            (
+                const_cast<Time&>(runTime),
+                nAlphaSubCycles
+            );
+            !(++alphaSubCycle).end();
+        )
+        {
+            solveAlphas();
+
+            int phasei = 0;
+            forAllIter(PtrDictionary<phaseModel>, phases(), iter)
+            {
+                phiSums[phasei] += iter().phi();
+                phasei++;
+            }
+        }
+
+        phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, phases(), iter)
+        {
+            phaseModel& phase = iter();
+            volScalarField& alpha = phase;
+
+            phase.phi() = phiSums[phasei]/nAlphaSubCycles;
+
+            // Correct the time index of the field
+            // to correspond to the global time
+            alpha.timeIndex() = runTime.timeIndex();
+
+            // Reset the old-time field value
+            alpha.oldTime() = alpha0s[phasei];
+            alpha.oldTime().timeIndex() = runTime.timeIndex();
+
+            phasei++;
+        }
+    }
+    else
+    {
+        solveAlphas();
+    }
+
+    forAllIter(PtrDictionary<phaseModel>, phases(), iter)
+    {
+        phaseModel& phase = iter();
+        phase.alphaRhoPhi() = fvc::interpolate(phase.rho())*phase.alphaPhi();
+    }
+
+    calcAlphas();
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2015 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::multiphaseSystem
+
+Description
+    Class which solves the volume fraction equations for two phases.
+
+SourceFiles
+    multiphaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseSystem_H
+#define multiphaseSystem_H
+
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class dragModel;
+class virtualMassModel;
+
+/*---------------------------------------------------------------------------*\
+                         Class multiphaseSystem Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseSystem
+:
+    public phaseSystem
+{
+    // Private data
+
+        volScalarField alphas_;
+
+        typedef HashTable<scalar, phasePairKey, phasePairKey::hash>
+            cAlphaTable;
+
+        cAlphaTable cAlphas_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        //- Conversion factor for degrees into radians
+        static const scalar convertToRad;
+
+
+    // Private member functions
+
+        void calcAlphas();
+
+        void solveAlphas();
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        void correctContactAngle
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2,
+            surfaceVectorField::GeometricBoundaryField& nHatb
+        ) const;
+
+        tmp<volScalarField> K
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2
+        ) const;
+
+
+        //- Return the drag coefficient for phase pair
+        virtual tmp<volScalarField> Kd(const phasePairKey& key) const = 0;
+
+        //- Return the face drag coefficient for phase pair
+        virtual tmp<surfaceScalarField> Kdf(const phasePairKey& key) const = 0;
+
+        //- Return the virtual mass coefficient for phase pair
+        virtual tmp<volScalarField> Vm(const phasePairKey& key) const = 0;
+
+        //- Return the face virtual mass coefficient for phase pair
+        virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const = 0;
+
+        //- Return the turbulent diffusivity for phase pair
+        //  Multiplies the phase-fraction gradient
+        virtual tmp<volScalarField> D(const phasePairKey& key) const = 0;
+
+        //- Return the interfacial mass flow rate for phase pair
+        virtual tmp<volScalarField> dmdt(const phasePairKey& key) const = 0;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseSystem");
+
+    // Declare runtime construction
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            multiphaseSystem,
+            dictionary,
+            (
+                const fvMesh& mesh
+            ),
+            (mesh)
+        );
+
+
+    // Constructors
+
+        //- Construct from fvMesh
+        multiphaseSystem(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~multiphaseSystem();
+
+
+    // Selectors
+
+        static autoPtr<multiphaseSystem> New
+        (
+            const fvMesh& mesh
+        );
+
+
+    // Member Functions
+
+        //- Return the drag coefficient for all phase-pairs
+        virtual const phaseSystem::KdTable& Kds() const = 0;
+
+        //- Return the drag coefficient for phase
+        virtual tmp<volScalarField> Kd(const phaseModel& phase) const = 0;
+
+        //- Return the combined force (lift + wall-lubrication) for phase pair
+        virtual autoPtr<PtrList<Foam::volVectorField> > Fs() const = 0;
+
+        //- Return the total interfacial mass transfer rate for phase
+        virtual tmp<volScalarField> dmdt(const phaseModel& phase) const = 0;
+
+        //- Return the momentum transfer matrices
+        virtual autoPtr<momentumTransferTable> momentumTransfer() const = 0;
+
+        //- Return the heat transfer matrices
+        virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
+
+        //- Return the mass transfer matrices
+        virtual autoPtr<massTransferTable> massTransfer() const = 0;
+
+        tmp<surfaceScalarField> surfaceTension(const phaseModel& phase) const;
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Solve for the phase fractions
+        virtual void solve();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "multiphaseSystemI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystemI.H

@@ -0,0 +1,28 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014-2015 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// ************************************************************************* //

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

@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 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 "addToRunTimeSelectionTable.H"
+
+#include "phaseSystem.H"
+#include "multiphaseSystem.H"
+#include "MomentumTransferPhaseSystem.H"
+#include "HeatTransferPhaseSystem.H"
+#include "InterfaceCompositionPhaseChangePhaseSystem.H"
+#include "ThermalPhaseChangePhaseSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef
+        HeatTransferPhaseSystem
+        <
+            MomentumTransferPhaseSystem<multiphaseSystem>
+        >
+        heatAndMomentumTransferMultiphaseSystem;
+
+    addNamedToRunTimeSelectionTable
+    (
+        multiphaseSystem,
+        heatAndMomentumTransferMultiphaseSystem,
+        dictionary,
+        heatAndMomentumTransferMultiphaseSystem
+    );
+
+    typedef
+        InterfaceCompositionPhaseChangePhaseSystem
+        <
+            MomentumTransferPhaseSystem<multiphaseSystem>
+        >
+        interfaceCompositionPhaseChangeMultiphaseSystem;
+
+    addNamedToRunTimeSelectionTable
+    (
+        multiphaseSystem,
+        interfaceCompositionPhaseChangeMultiphaseSystem,
+        dictionary,
+        interfaceCompositionPhaseChangeMultiphaseSystem
+    );
+
+    typedef
+        ThermalPhaseChangePhaseSystem
+        <
+            MomentumTransferPhaseSystem<multiphaseSystem>
+        >
+        thermalPhaseChangeMultiphaseSystem;
+
+    addNamedToRunTimeSelectionTable
+    (
+        multiphaseSystem,
+        thermalPhaseChangeMultiphaseSystem,
+        dictionary,
+        thermalPhaseChangeMultiphaseSystem
+    );
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 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 "multiphaseSystem.H"
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::multiphaseSystem> Foam::multiphaseSystem::New
+(
+    const fvMesh& mesh
+)
+{
+    const word multiphaseSystemType
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                propertiesName,
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).lookup("type")
+    );
+
+    Info<< "Selecting multiphaseSystem "
+        << multiphaseSystemType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(multiphaseSystemType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn("multiphaseSystem::New")
+            << "Unknown multiphaseSystemType type "
+            << multiphaseSystemType << endl << endl
+            << "Valid multiphaseSystem types are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return cstrIter()(mesh);
+}
+
+// ************************************************************************* //

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

@@ -0,0 +1,37 @@
+Info<< "Constructing momentum equations" << endl;
+
+PtrList<fvVectorMatrix> UEqns(fluid.phases().size());
+
+{
+    autoPtr<phaseSystem::momentumTransferTable>
+        momentumTransferPtr(fluid.momentumTransfer());
+
+    phaseSystem::momentumTransferTable&
+        momentumTransfer(momentumTransferPtr());
+
+    int phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        const volScalarField& alpha = phase;
+        const volScalarField& rho = phase.rho();
+        volVectorField& U = phase.U();
+
+        UEqns.set
+        (
+            phasei,
+            new fvVectorMatrix
+            (
+                phase.UEqn()
+             ==
+               *momentumTransfer[phase.name()]
+              + fvOptions(alpha, rho, U)
+            )
+        );
+
+        UEqns[phasei].relax();
+        fvOptions.constrain(UEqns[phasei]);
+
+        phasei++;
+    }
+}

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

@@ -0,0 +1,447 @@
+PtrList<surfaceScalarField> alphafs(fluid.phases().size());
+PtrList<volScalarField> rAUs(fluid.phases().size());
+PtrList<surfaceScalarField> alpharAUfs(fluid.phases().size());
+
+int phasei = 0;
+forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+{
+    phaseModel& phase = iter();
+    const volScalarField& alpha = phase;
+
+    alphafs.set(phasei, fvc::interpolate(alpha).ptr());
+    alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
+
+    rAUs.set
+    (
+        phasei,
+        new volScalarField
+        (
+            IOobject::groupName("rAU", phase.name()),
+            1.0
+           /(
+               UEqns[phasei].A()
+             + max(phase.residualAlpha() - alpha, scalar(0))
+              *phase.rho()/runTime.deltaT()
+            )
+        )
+    );
+
+    alpharAUfs.set
+    (
+        phasei,
+        (
+            fvc::interpolate(max(alpha, phase.residualAlpha())*rAUs[phasei])
+        ).ptr()
+    );
+
+    phasei++;
+}
+
+// Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
+PtrList<surfaceScalarField> phiFs(fluid.phases().size());
+{
+    autoPtr<PtrList<volVectorField> > Fs = fluid.Fs();
+
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+
+        phiFs.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                IOobject::groupName("phiF", phase.name()),
+                (fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
+            )
+        );
+
+        phasei++;
+    }
+}
+
+// --- Pressure corrector loop
+while (pimple.correct())
+{
+    // Update continuity errors due to temperature changes
+    fluid.correct();
+
+    PtrList<volVectorField> HbyAs(fluid.phases().size());
+
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        const volScalarField& alpha = phase;
+
+        // Correct fixed-flux BCs to be consistent with the velocity BCs
+        MRF.correctBoundaryFlux(phase.U(), phase.phi());
+
+        HbyAs.set
+        (
+            phasei,
+            new volVectorField
+            (
+                IOobject::groupName("HbyA", phase.name()),
+                phase.U()
+            )
+        );
+
+        HbyAs[phasei] =
+            rAUs[phasei]
+           *(
+                UEqns[phasei].H()
+              + max(phase.residualAlpha() - alpha, scalar(0))
+               *phase.rho()*phase.U().oldTime()/runTime.deltaT()
+            );
+
+        phasei++;
+    }
+
+    // Mean density for buoyancy force and p_rgh -> p
+    volScalarField rho("rho", fluid.rho());
+
+    surfaceScalarField ghSnGradRho
+    (
+        "ghSnGradRho",
+        ghf*fvc::snGrad(rho)*mesh.magSf()
+    );
+
+    PtrList<surfaceScalarField> phigs(fluid.phases().size());
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+
+        phigs.set
+        (
+            phasei,
+            (
+                alpharAUfs[phasei]
+               *(
+                   ghSnGradRho
+                 - (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
+                 - fluid.surfaceTension(phase)*mesh.magSf()
+                )
+            ).ptr()
+        );
+
+        phasei++;
+    }
+
+    PtrList<surfaceScalarField> phiHbyAs(fluid.phases().size());
+
+    surfaceScalarField phiHbyA
+    (
+        IOobject
+        (
+            "phiHbyA",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
+    );
+
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        const volScalarField& alpha = phase;
+
+        // ddtPhiCorr filter -- only apply in pure(ish) phases
+        surfaceScalarField alphafBar
+        (
+            fvc::interpolate(fvc::average(alphafs[phasei]))
+        );
+        surfaceScalarField phiCorrCoeff(pos(alphafBar - 0.99));
+
+        forAll(mesh.boundary(), patchi)
+        {
+            // Set ddtPhiCorr to 0 on non-coupled boundaries
+            if
+            (
+                !mesh.boundary()[patchi].coupled()
+             || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
+            )
+            {
+                phiCorrCoeff.boundaryField()[patchi] = 0;
+            }
+        }
+
+        phiHbyAs.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                IOobject::groupName("phiHbyA", phase.name()),
+                (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
+              + phiCorrCoeff
+               *fvc::interpolate
+                (
+                    alpha.oldTime()*phase.rho()().oldTime()*rAUs[phasei]
+                )
+               *(
+                   MRF.absolute(phase.phi().oldTime())
+                 - (fvc::interpolate(phase.U().oldTime()) & mesh.Sf())
+                )/runTime.deltaT()
+              - phiFs[phasei]
+              - phigs[phasei]
+            )
+        );
+
+        forAllConstIter
+        (
+            phaseSystem::KdTable,
+            fluid.Kds(),
+            KdIter
+        )
+        {
+            const volScalarField& K(*KdIter());
+
+            const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
+
+            const phaseModel* phase1 = &pair.phase1();
+            const phaseModel* phase2 = &pair.phase2();
+
+            forAllConstIter(phasePair, pair, iter)
+            {
+                if (phase1 == &phase)
+                {
+                    phiHbyAs[phasei] +=
+                        fvc::interpolate(rAUs[phasei]*K)
+                       *MRF.absolute(phase2->phi());
+
+                    HbyAs[phasei] += rAUs[phasei]*K*phase2->U();
+                }
+
+                Swap(phase1, phase2);
+            }
+        }
+
+        phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
+
+        phasei++;
+    }
+
+    MRF.makeRelative(phiHbyA);
+
+    // Construct pressure "diffusivity"
+    surfaceScalarField rAUf
+    (
+        IOobject
+        (
+            "rAUf",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
+    );
+
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        rAUf += alphafs[phasei]*alpharAUfs[phasei];
+        phasei++;
+    }
+    rAUf = mag(rAUf);
+
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    {
+        surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
+        phib = 0;
+        phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        {
+            phaseModel& phase = iter();
+            phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField();
+            phasei++;
+        }
+
+        setSnGrad<fixedFluxPressureFvPatchScalarField>
+        (
+            p_rgh.boundaryField(),
+            (
+                phiHbyA.boundaryField() - phib
+            )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+        );
+    }
+
+    PtrList<fvScalarMatrix> pEqnComps(fluid.phases().size());
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+
+        if (phase.compressible())
+        {
+            const volScalarField& alpha = phase;
+            const volScalarField& rho = phase.rho();
+
+            if (pimple.transonic())
+            {
+                surfaceScalarField phid
+                (
+                    IOobject::groupName("phid", phase.name()),
+                    fvc::interpolate(phase.thermo().psi())*phase.phi()
+                );
+
+                pEqnComps.set
+                (
+                    phasei,
+                    (
+                        (
+                            phase.continuityError() - fluid.dmdt(phase)
+                          - fvc::Sp
+                            (
+                                fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
+                                rho
+                            )
+                        )/rho
+                      + (alpha/rho)*correction
+                        (
+                            phase.thermo().psi()*fvm::ddt(p_rgh)
+                          + fvm::div(phid, p_rgh)
+                          - fvm::Sp(fvc::div(phid), p_rgh)
+                        )
+                    ).ptr()
+                );
+
+                deleteDemandDrivenData
+                (
+                    pEqnComps[phasei].faceFluxCorrectionPtr()
+                );
+                pEqnComps[phasei].relax();
+            }
+            else
+            {
+                pEqnComps.set
+                (
+                    phasei,
+                    (
+                        (
+                            phase.continuityError() - fluid.dmdt(phase)
+                          - fvc::Sp
+                            (
+                                (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())),
+                                rho
+                            )
+                        )/rho
+                      + (alpha*phase.thermo().psi()/rho)
+                       *correction(fvm::ddt(p_rgh))
+                    ).ptr()
+                );
+            }
+        }
+
+        phasei++;
+    }
+
+    // Cache p prior to solve for density update
+    volScalarField p_rgh_0(p_rgh);
+
+    // Iterate over the pressure equation to correct for non-orthogonality
+    while (pimple.correctNonOrthogonal())
+    {
+        // Construct the transport part of the pressure equation
+        fvScalarMatrix pEqnIncomp
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rAUf, p_rgh)
+        );
+
+        {
+            fvScalarMatrix pEqn(pEqnIncomp);
+
+            phasei = 0;
+            forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+            {
+                phaseModel& phase = iter();
+
+                if (phase.compressible())
+                {
+                    pEqn += pEqnComps[phasei];
+                }
+
+                phasei++;
+            }
+
+            solve
+            (
+                pEqn,
+                mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            );
+        }
+
+        // Correct fluxes and velocities on last non-orthogonal iteration
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqnIncomp.flux();
+
+            surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
+
+            phasei = 0;
+            forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+            {
+                phaseModel& phase = iter();
+
+                phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp;
+
+                // Set the phase dilatation rates
+                if (phase.compressible())
+                {
+                    phase.divU(-pEqnComps[phasei] & p_rgh);
+                }
+
+                phasei++;
+            }
+
+            // Optionally relax pressure for velocity correction
+            p_rgh.relax();
+
+            mSfGradp = pEqnIncomp.flux()/rAUf;
+
+            phasei = 0;
+            forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+            {
+                phaseModel& phase = iter();
+
+                phase.U() =
+                    HbyAs[phasei]
+                  + fvc::reconstruct
+                    (
+                        alpharAUfs[phasei]*mSfGradp
+                      - phiFs[phasei]
+                      - phigs[phasei]
+                    );
+                phase.U().correctBoundaryConditions();
+                fvOptions.correct(phase.U());
+
+                phasei++;
+            }
+        }
+    }
+
+    // Update and limit the static pressure
+    p = max(p_rgh + rho*gh, pMin);
+
+    // Limit p_rgh
+    p_rgh = p - rho*gh;
+
+    // Update densities from change in p_rgh
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        phase.rho()() += phase.thermo().psi()*(p_rgh - p_rgh_0);
+    }
+
+    // Correct p_rgh for consistency with p and the updated densities
+    rho = fluid.rho();
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+}

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

@@ -0,0 +1,285 @@
+// --- Pressure corrector loop
+while (pimple.correct())
+{
+    // rho1 = rho10 + psi1*p_rgh;
+    // rho2 = rho20 + psi2*p_rgh;
+
+    // tmp<fvScalarMatrix> pEqnComp1;
+    // tmp<fvScalarMatrix> pEqnComp2;
+
+    // //if (transonic)
+    // //{
+    // //}
+    // //else
+    // {
+    //     surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1);
+    //     surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
+
+    //     pEqnComp1 =
+    //         fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+    //       + fvc::div(phid1, p_rgh)
+    //       - fvc::Sp(fvc::div(phid1), p_rgh);
+
+    //     pEqnComp2 =
+    //         fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+    //       + fvc::div(phid2, p_rgh)
+    //       - fvc::Sp(fvc::div(phid2), p_rgh);
+    // }
+
+    PtrList<surfaceScalarField> alphafs(fluid.phases().size());
+    PtrList<volVectorField> HbyAs(fluid.phases().size());
+    PtrList<surfaceScalarField> phiHbyAs(fluid.phases().size());
+    PtrList<volScalarField> rAUs(fluid.phases().size());
+    PtrList<surfaceScalarField> rAlphaAUfs(fluid.phases().size());
+
+    int phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+
+        MRF.makeAbsolute(phase.phi().oldTime());
+        MRF.makeAbsolute(phase.phi());
+
+        HbyAs.set(phasei, new volVectorField(phase.U()));
+        phiHbyAs.set(phasei, new surfaceScalarField(1.0*phase.phi()));
+
+        phasei++;
+    }
+
+    surfaceScalarField phiHbyA
+    (
+        IOobject
+        (
+            "phiHbyA",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
+    );
+
+    volScalarField rho("rho", fluid.rho());
+    surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        const volScalarField& alpha = phase;
+
+        alphafs.set(phasei, fvc::interpolate(alpha).ptr());
+        alphafs[phasei].rename("hmm" + alpha.name());
+
+        volScalarField dragCoeffi
+        (
+            IOobject
+            (
+                "dragCoeffi",
+                runTime.timeName(),
+                mesh
+            ),
+            fluid.Kd(phase),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        dragCoeffi.correctBoundaryConditions();
+
+        rAUs.set(phasei, (1.0/(UEqns[phasei].A() + dragCoeffi)).ptr());
+        rAlphaAUfs.set
+        (
+            phasei,
+            (
+                alphafs[phasei]*fvc::interpolate(phase.rho())
+               /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
+            ).ptr()
+        );
+
+        HbyAs[phasei] = rAUs[phasei]*UEqns[phasei].H();
+
+        phiHbyAs[phasei] =
+        (
+            (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
+          + rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
+        );
+        MRF.makeRelative(phiHbyAs[phasei]);
+        MRF.makeRelative(phase.phi().oldTime());
+        MRF.makeRelative(phase.phi());
+
+        phiHbyAs[phasei] +=
+            rAlphaAUfs[phasei]
+           *(
+               fluid.surfaceTension(phase)*mesh.magSf()
+             + (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
+             - ghSnGradRho
+            )/fvc::interpolate(phase.rho());
+
+        forAllConstIter
+        (
+            phaseSystem::KdTable,
+            fluid.Kds(),
+            KdIter
+        )
+        {
+            const volScalarField& K(*KdIter());
+
+            const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
+
+            const phaseModel* phase1 = &pair.phase1();
+            const phaseModel* phase2 = &pair.phase2();
+
+            forAllConstIter(phasePair, pair, iter)
+            {
+                if (phase1 == &phase)
+                {
+                    phiHbyAs[phasei] +=
+                        fvc::interpolate(K)
+                       /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
+                       *phase2->phi();
+
+                    HbyAs[phasei] += rAUs[phasei]*K*phase2->U();
+                }
+
+                Swap(phase1, phase2);
+            }
+        }
+
+        phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
+
+        phasei++;
+    }
+
+    surfaceScalarField rAUf
+    (
+        IOobject
+        (
+            "rAUf",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
+    );
+
+    phasei = 0;
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        phaseModel& phase = iter();
+        rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])
+            /fvc::interpolate(phase.rho());
+
+        phasei++;
+    }
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    {
+        surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
+        phib = 0;
+        phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        {
+            phaseModel& phase = iter();
+
+            phib +=
+                alphafs[phasei].boundaryField()
+               *(mesh.Sf().boundaryField() & phase.U().boundaryField());
+
+            phasei++;
+        }
+
+        setSnGrad<fixedFluxPressureFvPatchScalarField>
+        (
+            p_rgh.boundaryField(),
+            (
+                phiHbyA.boundaryField() - MRF.relative(phib)
+            )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+        );
+    }
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqnIncomp
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rAUf, p_rgh)
+        );
+
+        pEqnIncomp.setReference(pRefCell, pRefValue);
+
+        solve
+        (
+            // (
+            //     (alpha1/rho1)*pEqnComp1()
+            //   + (alpha2/rho2)*pEqnComp2()
+            // ) +
+            pEqnIncomp,
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
+
+            phasei = 0;
+            phi = dimensionedScalar("phi", phi.dimensions(), 0);
+            forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+            {
+                phaseModel& phase = iter();
+
+                phase.phi() =
+                    phiHbyAs[phasei]
+                  + rAlphaAUfs[phasei]*mSfGradp/fvc::interpolate(phase.rho());
+                phi +=
+                    alphafs[phasei]*phiHbyAs[phasei]
+                  + mag(alphafs[phasei]*rAlphaAUfs[phasei])
+                   *mSfGradp/fvc::interpolate(phase.rho());
+
+                phasei++;
+            }
+
+            // dgdt =
+
+            // (
+            //     pos(alpha2)*(pEqnComp2 & p)/rho2
+            //   - pos(alpha1)*(pEqnComp1 & p)/rho1
+            // );
+
+            p_rgh.relax();
+
+            p = p_rgh + rho*gh;
+
+            mSfGradp = pEqnIncomp.flux()/rAUf;
+
+            phasei = 0;
+            forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+            {
+                phaseModel& phase = iter();
+
+                phase.U() =
+                    HbyAs[phasei]
+                  + fvc::reconstruct
+                    (
+                        rAlphaAUfs[phasei]
+                       *(
+                            fvc::interpolate(phase.rho() - rho)
+                           *(g & mesh.Sf())
+                          - ghSnGradRho
+                          + mSfGradp
+                        )
+                    )/phase.rho();
+
+                phase.U().correctBoundaryConditions();
+
+                phasei++;
+            }
+        }
+    }
+
+    //p = max(p, pMin);
+
+    // rho1 = rho10 + psi1*p_rgh;
+    // rho2 = rho20 + psi2*p_rgh;
+
+    // Dp1Dt = fvc::DDt(phi1, p_rgh);
+    // Dp2Dt = fvc::DDt(phi2, p_rgh);
+}

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

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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
+    reactingMultiphaseEulerFoam
+
+Description
+    Solver for a system of any number of compressible fluid phases with a
+    common pressure, but otherwise separate properties. The type of phase model
+    is run time selectable and can optionally represent multiple species and
+    in-phase reactions. The phase system is also run time selectable and can
+    optionally represent different types of momentun, heat and mass transfer.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseSystem.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+#include "pimpleControl.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
+    #include "createFields.H"
+
+    if (!LTS)
+    {
+        #include "CourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
+    // Switch faceMomentum
+    // (
+    //     pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
+    // );
+
+    // Switch implicitPhasePressure
+    // (
+    //     mesh.solverDict(alpha1.name()).lookupOrDefault<Switch>
+    //     (
+    //         "implicitPhasePressure", false
+    //     )
+    // );
+
+    //#include "pUf/createDDtU.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "CourantNo.H"
+            #include "setDeltaT.H"
+        }
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            fluid.solve();
+            fluid.correct();
+
+            #include "YEqns.H"
+
+            // if (faceMomentum)
+            // {
+            //     #include "pUf/UEqns.H"
+            //     #include "EEqns.H"
+            //     #include "pUf/pEqn.H"
+            //     #include "pUf/DDtU.H"
+            // }
+            // else
+            {
+                #include "pU/UEqns.H"
+                #include "EEqns.H"
+                #include "pU/pEqn.H"
+            }
+
+            fluid.correctKinematics();
+
+            if (pimple.turbCorr())
+            {
+                fluid.correctTurbulence();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = "
+            << runTime.elapsedCpuTime()
+            << " s\n\n" << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,36 @@
+{
+    volScalarField& rDeltaT = trDeltaT();
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        runTime.controlDict().lookupOrDefault<scalar>
+        (
+            "rDeltaTSmoothingCoeff",
+            0.02
+        )
+    );
+
+    surfaceScalarField maxPhi("maxPhi", phi);
+
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        maxPhi = max(maxPhi, mag(iter().phi()));
+    }
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.dimensionedInternalField() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(maxPhi)().dimensionedInternalField()
+       /((2*maxCo)*mesh.V())
+    );
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+}

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

@@ -50,8 +50,14 @@
 
 fluid.correctThermo();
 
-Info<< "  phase1.thermo().T(): " << min(phase1.thermo().T()).value()
+Info<< phase1.name() << " min/max T "
-    << " - " << max(phase1.thermo().T()).value() << endl;
+    << min(phase1.thermo().T()).value()
+    << " - "
+    << max(phase1.thermo().T()).value()
+    << endl;
 
-Info<< "  phase2.thermo().T(): " << min(phase2.thermo().T()).value()
+Info<< phase2.name() << " min/max T "
-    << " - " << max(phase2.thermo().T()).value() << endl;
+    << min(phase2.thermo().T()).value()
+    << " - "
+    << max(phase2.thermo().T()).value()
+    << endl;

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

@@ -237,58 +237,72 @@ while (pimple.correct())
     // Construct the compressibility parts of the pressure equation
     if (pimple.transonic())
     {
-        surfaceScalarField phid1
+        if (phase1.compressible())
-        (
+        {
-            IOobject::groupName("phid", phase1.name()),
+            surfaceScalarField phid1
-            fvc::interpolate(psi1)*phi1
-        );
-        surfaceScalarField phid2
-        (
-            IOobject::groupName("phid", phase2.name()),
-            fvc::interpolate(psi2)*phi2
-        );
-
-        pEqnComp1 =
             (
-                phase1.continuityError() - fluid.dmdt()
+                IOobject::groupName("phid", phase1.name()),
-              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+                fvc::interpolate(psi1)*phi1
-            )/rho1
-          + (alpha1/rho1)*correction
-            (
-                psi1*fvm::ddt(p_rgh)
-              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
             );
-        deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
-        pEqnComp1().relax();
 
-        pEqnComp2 =
+            pEqnComp1 =
+                (
+                    phase1.continuityError() - fluid.dmdt()
+                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+                )/rho1
+              + (alpha1/rho1)*correction
+                (
+                    psi1*fvm::ddt(p_rgh)
+                  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+                );
+
+            deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
+            pEqnComp1().relax();
+        }
+
+        if (phase2.compressible())
+        {
+            surfaceScalarField phid2
             (
-                phase2.continuityError() + fluid.dmdt()
+                IOobject::groupName("phid", phase2.name()),
-              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+                fvc::interpolate(psi2)*phi2
-            )/rho2
-          + (alpha2/rho2)*correction
-            (
-                psi2*fvm::ddt(p_rgh)
-              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
             );
-        deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
+
-        pEqnComp2().relax();
+            pEqnComp2 =
+                (
+                    phase2.continuityError() + fluid.dmdt()
+                  - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+                )/rho2
+               + (alpha2/rho2)*correction
+                (
+                    psi2*fvm::ddt(p_rgh)
+                  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+                );
+            deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
+            pEqnComp2().relax();
+        }
     }
     else
     {
-        pEqnComp1 =
+        if (phase1.compressible())
-            (
+        {
-                phase1.continuityError() - fluid.dmdt()
+            pEqnComp1 =
-              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+                (
-            )/rho1
+                    phase1.continuityError() - fluid.dmdt()
-          + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
+                  - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+                )/rho1
+              + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
+        }
 
-        pEqnComp2 =
+        if (phase2.compressible())
-            (
+        {
-                phase2.continuityError() + fluid.dmdt()
+            pEqnComp2 =
-              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+                (
-            )/rho2
+                    phase2.continuityError() + fluid.dmdt()
-          + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
+                  - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+                )/rho2
+              + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
+        }
     }
 
     // Cache p prior to solve for density update
@@ -304,11 +318,25 @@ while (pimple.correct())
           - fvm::laplacian(rAUf, p_rgh)
         );
 
-        solve
+        {
-        (
+            fvScalarMatrix pEqn(pEqnIncomp);
-            pEqnComp1() + pEqnComp2() + pEqnIncomp,
+
-            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            if (phase1.compressible())
-        );
+            {
+                pEqn += pEqnComp1();
+            }
+
+            if (phase2.compressible())
+            {
+                pEqn += pEqnComp2();
+            }
+
+            solve
+            (
+                pEqn,
+                mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+            );
+        }
 
         // Correct fluxes and velocities on last non-orthogonal iteration
         if (pimple.finalNonOrthogonalIter())

tutorials/multiphase/multiphaseEulerFoam/bubbleColumn/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     "div\(phi,alpha.*\)"    Gauss vanLeer;
     "div\(phir,alpha.*,alpha.*\)"   Gauss vanLeer;
 
-    "div\(phiAlpha.*,U.*\)" Gauss limitedLinearV 1;
+    "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
     div(Rc)                 Gauss linear;
     "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
 }

tutorials/multiphase/multiphaseEulerFoam/damBreak4phase/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
     "div\(phi,alpha.*\)"    Gauss vanLeer;
     "div\(phir,alpha.*,alpha.*\)"   Gauss vanLeer;
 
-    "div\(phiAlpha.*,U.*\)" Gauss limitedLinearV 1;
+    "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
     div(Rc)                 Gauss linear;
     "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
 }

tutorials/multiphase/multiphaseEulerFoam/damBreak4phaseFine/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
     "div\(phi,alpha.*\)"    Gauss vanLeer;
     "div\(phir,alpha.*,alpha.*\)"   Gauss vanLeer;
 
-    "div\(phiAlpha.*,U.*\)" Gauss limitedLinearV 1;
+    "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
     div(Rc)                 Gauss linear;
     "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
 }

tutorials/multiphase/multiphaseEulerFoam/mixerVessel2D/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
     "div\(phi,alpha.*\)"    Gauss vanLeer;
     "div\(phir,alpha.*,alpha.*\)"   Gauss vanLeer;
 
-    "div\(phiAlpha.*,U.*\)" Gauss limitedLinearV 1;
+    "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
     div(Rc)                 Gauss linear;
     "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
 }

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/T.air

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Tair;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 300;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.air;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/T.water

@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Twater;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 0 0 1 0 0 0];
+
+internalField       uniform 350;
+
+boundaryField
+{
+    walls
+    {
+        type               zeroGradient;
+    }
+    outlet
+    {
+        type               inletOutlet;
+        phi                phi.water;
+        inletValue         $internalField;
+        value              $internalField;
+    }
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    frontAndBackPlanes
+    {
+        type               empty;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/Theta

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Theta;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 1.0e-7;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1.0e-7;
+        value           uniform 1.0e-7;
+    }
+
+    walls
+    {
+        type            zeroGradient;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/U.air

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0.1 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.air;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/U.water

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedValue;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               pressureInletOutletVelocity;
+        phi                phi.water;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedValue;
+        value              uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/alpha.air.org

@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.5;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.air;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/alpha.water.org

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.5;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        phi             phi.water;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    walls
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/p

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               calculated;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               calculated;
+        value              $internalField;
+    }
+    walls
+    {
+        type               calculated;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/p_rgh

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions          [1 -1 -2 0 0 0 0];
+
+internalField       uniform 1e5;
+
+boundaryField
+{
+    inlet
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+    outlet
+    {
+        type               prghPressure;
+        p                  $internalField;
+        value              $internalField;
+    }
+    walls
+    {
+        type               fixedFluxPressure;
+        value              $internalField;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties

@@ -0,0 +1,169 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      phaseProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+type heatAndMomentumTransferMultiphaseSystem;
+
+phases (air water);
+
+air
+{
+    type            purePhaseModel;
+    diameterModel   isothermal;
+    isothermalCoeffs
+    {
+        d0              3e-3;
+        p0              1e5;
+    }
+
+    residualAlpha   1e-6;
+}
+
+water
+{
+    type            purePhaseModel;
+    diameterModel   constant;
+    constantCoeffs
+    {
+        d               1e-4;
+    }
+
+    residualAlpha   1e-6;
+}
+
+blending
+{
+    default
+    {
+        type            linear;
+        minFullyContinuousAlpha.air 0.7;
+        minPartlyContinuousAlpha.air 0.5;
+        minFullyContinuousAlpha.water 0.7;
+        minPartlyContinuousAlpha.water 0.5;
+    }
+}
+
+surfaceTension
+(
+    (air and water)
+    {
+        type            constant;
+        sigma           0.07;
+    }
+);
+
+interfaceCompression
+();
+
+aspectRatio
+(
+    (air in water)
+    {
+        type            constant;
+        E0              1.0;
+    }
+
+    (water in air)
+    {
+        type            constant;
+        E0              1.0;
+    }
+);
+
+drag
+(
+    (air in water)
+    {
+        type            SchillerNaumann;
+
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (water in air)
+    {
+        type            SchillerNaumann;
+
+        residualRe      1e-3;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+
+    (air and water)
+    {
+        type            segregated;
+
+        m               0.5;
+        n               8;
+        swarmCorrection
+        {
+            type        none;
+        }
+    }
+);
+
+virtualMass
+(
+    (air in water)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+
+    (water in air)
+    {
+        type            constantCoefficient;
+        Cvm             0.5;
+    }
+);
+
+heatTransfer
+(
+    (air in water)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+
+    (water in air)
+    {
+        type            RanzMarshall;
+        residualAlpha   1e-4;
+    }
+);
+
+lift
+(
+);
+
+wallLubrication
+(
+);
+
+turbulentDispersion
+(
+);
+
+// Minimum allowable pressure
+pMin            10000;
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/polyMesh/boundary

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+4
+(
+    inlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       3650;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          25;
+        startFace       3675;
+    }
+    walls
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          150;
+        startFace       3700;
+    }
+    defaultFaces
+    {
+        type            empty;
+        inGroups        1(empty);
+        nFaces          3750;
+        startFace       3850;
+    }
+)
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/thermophysicalProperties.air

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cp          1007;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.84e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/thermophysicalProperties.water

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectFluid;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles      1;
+        molWeight   18;
+    }
+    equationOfState
+    {
+        R           3000;
+        rho0        1027;
+    }
+    thermodynamics
+    {
+        Cp          4195;
+        Hf          0;
+    }
+    transport
+    {
+        mu          3.645e-4;
+        Pr          2.289;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/turbulenceProperties.air

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/turbulenceProperties.water

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/system/blockMeshDict

@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 0)
+    (0.15 0 0)
+    (0.15 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (0.15 0 0.1)
+    (0.15 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (25 75 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    patch inlet
+    (
+        (1 5 4 0)
+    )
+    patch outlet
+    (
+        (3 7 6 2)
+    )
+    wall walls
+    (
+        (0 4 7 3)
+        (2 6 5 1)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/system/controlDict

@@ -0,0 +1,95 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     reactingMultiphaseEulerFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          0.005;
+
+writeControl    runTime;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  no;
+
+maxCo           0.5;
+
+maxDeltaT       1;
+
+functions0
+{
+    fieldAverage1
+    {
+        type            fieldAverage;
+        functionObjectLibs ( "libfieldFunctionObjects.so" );
+        outputControl   outputTime;
+        fields
+        (
+            U.air
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            U.water
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            alpha.air
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+
+            p
+            {
+                 mean        on;
+                 prime2Mean  off;
+                 base        time;
+            }
+        );
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/system/fvSchemes

@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default                         none;
+
+    div(phi,alpha.air)              Gauss vanLeer;
+    div(phi,alpha.water)            Gauss vanLeer;
+    div(phir,alpha.water,alpha.air) Gauss vanLeer;
+    div(phir,alpha.air,alpha.water) Gauss vanLeer;
+
+    "div\(alphaRhoPhi.*,U.*\)"      Gauss limitedLinearV 1;
+    "div\(phi.*,U.*\)"              Gauss limitedLinearV 1;
+
+    "div\(alphaRhoPhi.*,(h|e).*\)"  Gauss limitedLinear 1;
+    "div\(alphaRhoPhi.*,K.*\)"      Gauss limitedLinear 1;
+    "div\(alphaPhi.*,p\)"           Gauss limitedLinear 1;
+
+    "div\(\(\(\(alpha.*\*thermo:rho.*\)*nuEff.*\)\*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear uncorrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         uncorrected;
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/system/fvSolution

@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "alpha.*"
+    {
+        nAlphaCorr      1;
+        nAlphaSubCycles 2;
+    }
+
+    p_rgh
+    {
+        solver          GAMG;
+        smoother        DIC;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        nFinestSweeps   2;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+        tolerance       1e-8;
+        relTol          0;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    "U.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-5;
+        relTol          0;
+        minIter         1;
+    }
+
+    "e.*"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
+        minIter         1;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 3;
+    nCorrectors      1;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/system/setFieldsDict

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha.air 1
+);
+
+regions
+(
+    boxToCell
+    {
+        box (0 0 -0.1) (0.15 0.701 0.1);
+        fieldValues
+        (
+            volScalarFieldValue alpha.air 0
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/U

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/U.air

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/U.mercury

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U.mercury;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/U.oil

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U.oil;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/U.water

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      Uwater;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    leftWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    rightWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    lowerWall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    atmosphere
+    {
+        type            fluxCorrectedVelocity;
+        value           uniform (0 0 0);
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/alpha.air

@@ -0,0 +1,79 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 0;
+    }
+    rightWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 1;
+    }
+    lowerWall
+    {
+        type            alphaContactAngle;
+        thetaProperties
+        (
+            ( water air ) 90 0 0 0
+            ( oil air ) 90 0 0 0
+            ( mercury air ) 90 0 0 0
+            ( water oil ) 90 0 0 0
+            ( water mercury ) 90 0 0 0
+            ( oil mercury ) 90 0 0 0
+        );
+        value           uniform 0;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/alpha.mercury

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.mercury;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/alpha.oil

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.oil;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/alpha.water

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            inletOutlet;
+        inletValue      uniform 0;
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/damBreak4phase/0.org/alphas

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    leftWall
+    {
+        type            zeroGradient;
+    }
+    rightWall
+    {
+        type            zeroGradient;
+    }
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+    atmosphere
+    {
+        type            zeroGradient;
+    }
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //