Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Conflicts:
	tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/p
	tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/changeDictionaryDict
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/H2O
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/T
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/U
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/air
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/alphat
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/k
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/mut
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/omega
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/0/p
	tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/polyMesh/boundary

applications/solvers/combustion/PDRFoam/EaEqn.H

@@ -0,0 +1,17 @@
+{
+    volScalarField& hea = thermo.he();
+
+    solve
+    (
+        betav*fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea)
+      + betav*fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            hea.name() == "ea"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -betav*dpdt
+        )
+      - fvm::laplacian(Db, hea)
+    );
+
+    thermo.correct();
+}

applications/solvers/combustion/PDRFoam/EauEqn.H

@@ -0,0 +1,22 @@
+if (ign.ignited())
+{
+    volScalarField& heau = thermo.heu();
+
+    solve
+    (
+        betav*fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
+      + (betav*fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
+      + (
+            heau.name() == "eau"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))*rho/thermo.rhou()
+          : -betav*dpdt*rho/thermo.rhou()
+        )
+      - fvm::laplacian(Db, heau)
+
+        // These terms cannot be used in partially-premixed combustion due to
+        // the resultant inconsistency between ft and heau transport.
+        // A possible solution would be to solve for ftu as well as ft.
+        //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
+        //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
+    );
+}

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -123,12 +123,12 @@ int main(int argc, char *argv[])
             {
                 #include "bEqn.H"
                 #include "ftEqn.H"
-                #include "hauEqn.H"
+                #include "EauEqn.H"
-                #include "haEqn.H"
+                #include "EaEqn.H"
 
                 if (!ign.ignited())
                 {
-                    hau == ha;
+                    thermo.heu() == thermo.he();
                 }
 
                 #include "pEqn.H"

applications/solvers/combustion/PDRFoam/createFields.H

@@ -5,6 +5,8 @@
         psiuReactionThermo::New(mesh)
     );
     psiuReactionThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "ha", "ea");
+
     basicMultiComponentMixture& composition = thermo.composition();
 
     volScalarField rho
@@ -22,14 +24,10 @@
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
-    volScalarField& ha = thermo.he();
-    volScalarField& hau = thermo.heu();
 
     volScalarField& b = composition.Y("b");
     Info<< "min(b) = " << min(b).value() << endl;
 
-    //const volScalarField& T = thermo->T();
-
     Info<< "\nReading field U\n" << endl;
     volVectorField U
     (
@@ -196,6 +194,6 @@
     }
 
     fields.add(b);
-    fields.add(ha);
+    fields.add(thermo.he());
-    fields.add(hau);
+    fields.add(thermo.heu());
     flameWrinkling->addXi(fields);

applications/solvers/combustion/PDRFoam/haEqn.H

@@ -1,13 +0,0 @@
-{
-    solve
-    (
-        betav*fvm::ddt(rho, ha)
-      + mvConvection->fvmDiv(phi, ha)
-      - fvm::laplacian(Db, ha)
-     ==
-        betav*dpdt
-      - betav*(fvc::ddt(rho, K) + fvc::div(phi, K))
-    );
-
-    thermo.correct();
-}

applications/solvers/combustion/PDRFoam/hauEqn.H

@@ -1,18 +0,0 @@
-if (ign.ignited())
-{
-    solve
-    (
-        betav*fvm::ddt(rho, hau)
-      + mvConvection->fvmDiv(phi, hau)
-      - fvm::laplacian(Db, hau)
-
-    // These terms cannot be used in partially-premixed combustion due to
-    // the resultant inconsistency between ft and hau transport.
-    // A possible solution would be to solve for ftu as well as ft.
-    //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), hau)
-    //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hau)
-
-     ==
-        betav*(dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou()
-    );
-}

applications/solvers/combustion/XiFoam/EaEqn.H

@@ -0,0 +1,20 @@
+{
+    volScalarField& hea = thermo.he();
+
+    fvScalarMatrix EaEqn
+    (
+        fvm::ddt(rho, hea) + mvConvection->fvmDiv(phi, hea)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            hea.name() == "ea"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), hea)
+    );
+
+    EaEqn.relax();
+    EaEqn.solve();
+
+    thermo.correct();
+}

applications/solvers/combustion/XiFoam/EauEqn.H

@@ -0,0 +1,22 @@
+if (ign.ignited())
+{
+    volScalarField& heau = thermo.heu();
+
+    solve
+    (
+        fvm::ddt(rho, heau) + mvConvection->fvmDiv(phi, heau)
+      + (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
+      + (
+            heau.name() == "eau"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))*rho/thermo.rhou()
+          : -dpdt*rho/thermo.rhou()
+        )
+      - fvm::laplacian(turbulence->alphaEff(), heau)
+
+        // These terms cannot be used in partially-premixed combustion due to
+        // the resultant inconsistency between ft and heau transport.
+        // A possible solution would be to solve for ftu as well as ft.
+        //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
+        //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
+    );
+}

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -97,12 +97,12 @@ int main(int argc, char *argv[])
 
             #include "ftEqn.H"
             #include "bEqn.H"
-            #include "hauEqn.H"
+            #include "EauEqn.H"
-            #include "haEqn.H"
+            #include "EaEqn.H"
 
             if (!ign.ignited())
             {
-                hau == ha;
+                thermo.heu() == thermo.he();
             }
 
             // --- Pressure corrector loop

applications/solvers/combustion/XiFoam/createFields.H

@@ -5,6 +5,8 @@
         psiuReactionThermo::New(mesh)
     );
     psiuReactionThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "ha", "ea");
+
     basicMultiComponentMixture& composition = thermo.composition();
 
     volScalarField rho
@@ -22,14 +24,10 @@
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
-    volScalarField& ha = thermo.he();
-    volScalarField& hau = thermo.heu();
 
     volScalarField& b = composition.Y("b");
     Info<< "min(b) = " << min(b).value() << endl;
 
-    const volScalarField& T = thermo.T();
-
 
     Info<< "\nReading field U\n" << endl;
     volVectorField U
@@ -138,5 +136,5 @@
     }
 
     fields.add(b);
-    fields.add(ha);
+    fields.add(thermo.he());
-    fields.add(hau);
+    fields.add(thermo.heu());

applications/solvers/combustion/XiFoam/haEqn.H

@@ -1,16 +0,0 @@
-{
-    fvScalarMatrix haEqn
-    (
-        fvm::ddt(rho, ha)
-      + mvConvection->fvmDiv(phi, ha)
-      - fvm::laplacian(turbulence->alphaEff(), ha)
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-    );
-
-    haEqn.relax();
-    haEqn.solve();
-
-    thermo.correct();
-}

applications/solvers/combustion/XiFoam/hauEqn.H

@@ -1,18 +0,0 @@
-if (ign.ignited())
-{
-    solve
-    (
-        fvm::ddt(rho, hau)
-      + mvConvection->fvmDiv(phi, hau)
-      - fvm::laplacian(turbulence->alphaEff(), hau)
-
-    // These terms cannot be used in partially-premixed combustion due to
-    // the resultant inconsistency between ft and hau transport.
-    // A possible solution would be to solve for ftu as well as ft.
-    //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), hau)
-    //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hau)
-
-    ==
-        (dpdt - (fvc::ddt(rho, K) + fvc::div(phi, K)))*rho/thermo.rhou()
-    );
-}

applications/solvers/combustion/chemFoam/createFields.H

@@ -30,6 +30,8 @@
     scalar dtChem = refCast<const psiChemistryModel>(chemistry).deltaTChem()[0];
 
     psiReactionThermo& thermo = chemistry.thermo();
+    thermo.validate(args.executable(), "h");
+
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
 
@@ -47,7 +49,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.he();
 
     volScalarField Rspecific
     (

applications/solvers/combustion/chemFoam/hEqn.H

@@ -1,10 +1,12 @@
 {
+    volScalarField& h = thermo.he();
+
     if (constProp == "volume")
     {
-        hs[0] = u0 + p[0]/rho[0] + integratedHeat;
+        h[0] = u0 + p[0]/rho[0] + integratedHeat;
     }
     else
     {
-        hs[0] = hs0 + integratedHeat;
+        h[0] = h0 + integratedHeat;
     }
 }

applications/solvers/combustion/chemFoam/readInitialConditions.H

@@ -83,20 +83,19 @@
         }
     }
 
-    scalar hs0 = 0.0;
+    scalar h0 = 0.0;
     forAll(Y, i)
     {
         Y[i] = Y0[i];
-        hs0 += Y0[i]*specieData[i].Hs(p[i], T0);
+        h0 += Y0[i]*specieData[i].Hs(p[i], T0);
     }
 
-    hs = dimensionedScalar("h", dimEnergy/dimMass, hs0);
+    thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);
-
     thermo.correct();
 
     rho = thermo.rho();
     scalar rho0 = rho[0];
-    scalar u0 = hs0 - p0/rho0;
+    scalar u0 = h0 - p0/rho0;
     scalar R0 = p0/(rho0*T0);
     Rspecific[0] = R0;
 

applications/solvers/combustion/coldEngineFoam/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I../engineFoam \
     -I../XiFoam \
+    -I../../compressible/rhoPimpleFoam \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -81,7 +81,7 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                #include "hEqn.H"
+                #include "EEqn.H"
                 #include "pEqn.H"
             }
 

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -5,6 +5,7 @@
         psiThermo::New(mesh)
     );
     psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField rho
     (
@@ -21,7 +22,6 @@
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
-    volScalarField& h = thermo.he();
     const volScalarField& T = thermo.T();
 
 

applications/solvers/combustion/coldEngineFoam/hEqn.H

@@ -1,12 +0,0 @@
-{
-    solve
-    (
-        fvm::ddt(rho, h)
-      + fvm::div(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-      - fvc::div(phi, 0.5*magSqr(U))
-    );
-
-    thermo.correct();
-}

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -103,12 +103,12 @@ int main(int argc, char *argv[])
 
             #include "ftEqn.H"
             #include "bEqn.H"
-            #include "hauEqn.H"
+            #include "EauEqn.H"
-            #include "haEqn.H"
+            #include "EaEqn.H"
 
             if (!ign.ignited())
             {
-                hau == ha;
+                thermo.heu() == thermo.he();
             }
 
             // --- Pressure corrector loop

applications/solvers/combustion/engineFoam/logSummary.H

@@ -1,5 +1,6 @@
 Info<< "Mean pressure:" << p.weightedAverage(mesh.V()).value() << endl;
-Info<< "Mean temperature:" << T.weightedAverage(mesh.V()).value() << endl;
+Info<< "Mean temperature:" << thermo.T().weightedAverage(mesh.V()).value()
+    << endl;
 Info<< "Mean u':"
     << (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value()
     << endl;
@@ -7,7 +8,7 @@ Info<< "Mean u':"
 logSummaryFile
     << runTime.theta() << tab
     << p.weightedAverage(mesh.V()).value() << tab
-    << T.weightedAverage(mesh.V()).value() << tab
+    << thermo.T().weightedAverage(mesh.V()).value() << tab
     << (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value()
     << tab
     << 1 - b.weightedAverage(mesh.V()).value()

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -47,22 +47,27 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
 
-    fvScalarMatrix hsEqn
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
     (
-        fvm::ddt(rho, hs)
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
-      + mvConvection->fvmDiv(phi, hs)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        dpdt
+        combustion->Sh()
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + combustion->Sh()
       + radiation->Sh(thermo)
-      + parcels.Sh(hs)
+      + parcels.Sh(he)
       + surfaceFilm.Sh()
     );
 
-    hsEqn.relax();
+    EEqn.relax();
-    hsEqn.solve();
+    EEqn.solve();
 
     thermo.correct();
 

applications/solvers/combustion/fireFoam/createFields.H

@@ -11,6 +11,7 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
     psiReactionThermo& thermo = combustion->thermo();
+    thermo.validate(args.executable(), "h", "e");
 
     SLGThermo slgThermo(mesh, thermo);
 
@@ -34,7 +35,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.he();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -128,7 +128,7 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(hs);
+    fields.add(thermo.he());
 
     IOdictionary additionalControlsDict
     (

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -94,7 +94,7 @@ int main(int argc, char *argv[])
             while (pimple.loop())
             {
                 #include "UEqn.H"
-                #include "YhsEqn.H"
+                #include "YEEqn.H"
 
                 // --- Pressure corrector loop
                 while (pimple.correct())

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -0,0 +1,26 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+//    - fvm::laplacian(turbulence->muEff(), he)  // unit lewis no.
+     ==
+        reaction->Sh()
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+
+    Info<< "min/max(T) = "
+        << min(T).value() << ", " << max(T).value() << endl;
+}

applications/solvers/combustion/reactingFoam/createFields.H

@@ -6,6 +6,7 @@ autoPtr<combustionModels::psiCombustionModel> reaction
 );
 
 psiReactionThermo& thermo = reaction->thermo();
+thermo.validate(args.executable(), "h", "e");
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,7 +41,6 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& hs = thermo.he();
 const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
@@ -84,7 +84,7 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(hs);
+fields.add(thermo.he());
 
 volScalarField dQ
 (

applications/solvers/combustion/reactingFoam/hsEqn.H

@@ -1,21 +0,0 @@
-{
-    fvScalarMatrix hsEqn
-    (
-        fvm::ddt(rho, hs)
-      + mvConvection->fvmDiv(phi, hs)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
-//      - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + reaction->Sh()
-    );
-
-    hsEqn.relax();
-    hsEqn.solve();
-
-    thermo.correct();
-
-    Info<< "T gas min/max   = " << min(T).value() << ", "
-        << max(T).value() << endl;
-}

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,7 +70,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hsEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/combustion/rhoReactingFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../reactingFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \

applications/solvers/combustion/rhoReactingFoam/UEqn.H

@@ -1,16 +0,0 @@
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(phi, U)
-      + turbulence->divDevRhoReff(U)
-     ==
-        rho*g
-    );
-
-    UEqn.relax();
-
-    if (pimple.momentumPredictor())
-    {
-        solve(UEqn == -fvc::grad(p));
-        K = 0.5*magSqr(U);
-    }

applications/solvers/combustion/rhoReactingFoam/YEqn.H

@@ -1,47 +0,0 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
-(
-    fv::convectionScheme<scalar>::New
-    (
-        mesh,
-        fields,
-        phi,
-        mesh.divScheme("div(phi,Yi_h)")
-    )
-);
-
-{
-    reaction->correct();
-    dQ = reaction->dQ();
-    label inertIndex = -1;
-    volScalarField Yt(0.0*Y[0]);
-
-    forAll(Y, i)
-    {
-        if (Y[i].name() != inertSpecie)
-        {
-            volScalarField& Yi = Y[i];
-
-            fvScalarMatrix YiEqn
-            (
-                fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
-             ==
-                reaction->R(Yi)
-            );
-
-            YiEqn.relax();
-            YiEqn.solve(mesh.solver("Yi"));
-
-            Yi.max(0.0);
-            Yt += Yi;
-        }
-        else
-        {
-            inertIndex = i;
-        }
-    }
-
-    Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
-}

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -6,6 +6,7 @@ autoPtr<combustionModels::rhoCombustionModel> reaction
 );
 
 rhoReactionThermo& thermo = reaction->thermo();
+thermo.validate(args.executable(), "h", "e");
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,7 +41,6 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& hs = thermo.he();
 const volScalarField& T = thermo.T();
 
 
@@ -86,7 +86,7 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(hs);
+fields.add(thermo.he());
 
 volScalarField dQ
 (

applications/solvers/combustion/rhoReactingFoam/hsEqn.H

@@ -1,21 +0,0 @@
-{
-    fvScalarMatrix hsEqn
-    (
-        fvm::ddt(rho, hs)
-      + mvConvection->fvmDiv(phi, hs)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
-//    - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + reaction->Sh()
-    );
-
-    hsEqn.relax();
-    hsEqn.solve();
-
-    thermo.correct();
-
-   Info<< "min/max(T) = "
-        << min(T).value() << ", " << max(T).value() << endl;
-}

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -72,7 +72,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hsEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -4,7 +4,6 @@ tmp<fvVectorMatrix> UEqn
 (
     fvm::ddt(rho, U)
   + fvm::div(phi, U)
-  - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), U)
   + turbulence->divDevRhoReff(U)
 );
 

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -5,9 +5,9 @@
         psiThermo::New(mesh)
     );
     psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.he();
     const volScalarField& psi = thermo.psi();
 
     volScalarField rho

applications/solvers/compressible/rhoPimpleFoam/hEqn.H

@@ -1,20 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + fvm::div(phi, h)
-      - fvm::Sp(fvc::ddt(rho) + fvc::div(phi), h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        dpdt
-      - (
-            fvc::ddt(rho, K) + fvc::div(phi, K)
-          - (fvc::ddt(rho) + fvc::div(phi))*K
-        )
-    );
-
-    hEqn.relax();
-    hEqn.solve();
-
-    thermo.correct();
-}

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/rhoPimplecFoam.C

@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C

@@ -85,7 +85,7 @@ int main(int argc, char *argv[])
             turbulence->correct();
 
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C

@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "UEqn.H"
-            #include "hEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/compressible/rhoSimpleFoam/UEqn.H

@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H

@@ -7,7 +7,6 @@
     fvScalarMatrix TEqn
     (
         fvm::div(phi, T)
-      - fvm::Sp(fvc::div(phi), T)
       - fvm::laplacian(kappaEff, T)
     );
 

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../../compressible/rhoPimpleFoam \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/EEqn.H

@@ -0,0 +1,26 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turb.alphaEff(), he)
+     ==
+        rad.Sh(thermo)
+      + sources(rho, he)
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+    rad.correct();
+
+    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+        << max(thermo.T()).value() << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H

@@ -1,23 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::div(phi, h)
-      - fvm::Sp(fvc::div(phi), h)
-      - fvm::laplacian(turb.alphaEff(), h)
-     ==
-      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
-      + rad.Sh(thermo)
-      + sources(rho, h)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    rad.correct();
-
-    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-        << max(thermo.T()).value() << endl;
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -7,23 +7,31 @@
     volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    bool closedVolume = adjustPhi(phi, U, p_rgh);
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
+
+    phiHbyA += phig;
+
     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
     bool compressible = (compressibility.value() > SMALL);
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    phi += phig;
-
     // Solve pressure
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference
@@ -37,14 +45,14 @@
         if (nonOrth == nNonOrthCorr)
         {
             // Calculate the conservative fluxes
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -1,6 +1,8 @@
     const fvMesh& mesh = fluidRegions[i];
 
     rhoThermo& thermo = thermoFluid[i];
+    thermo.validate(args.executable(), "h", "e");
+
     volScalarField& rho = rhoFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
@@ -9,7 +11,6 @@
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
-    volScalarField& h = thermo.he();
 
     IObasicSourceList& sources = heatSources[i];
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H

@@ -4,7 +4,7 @@
     rho.storePrevIter();
     {
         #include "UEqn.H"
-        #include "hEqn.H"
+        #include "EEqn.H"
         #include "pEqn.H"
     }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H

@@ -2,7 +2,6 @@
     fvScalarMatrix hPorousEqn
     (
         fvm::div(porousPhi, porousH)
-      - fvm::Sp(fvc::div(porousPhi), porousH)
       - fvm::laplacian(turbPorous.alphaEff(), porousH)
      ==
       - fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)")

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/EEqn.H

@@ -0,0 +1,27 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + fvm::div(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turb.alphaEff(), he)
+     ==
+        rad.Sh(thermo)
+      + sources(rho, he)
+    );
+
+    EEqn.relax();
+    EEqn.solve(mesh.solver(he.select(finalIter)));
+
+    thermo.correct();
+    rad.correct();
+
+    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+        << max(thermo.T()).value() << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H

@@ -1,24 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + fvm::div(phi, h)
-      - fvm::laplacian(turb.alphaEff(), h)
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + rad.Sh(thermo)
-      + sources(rho, h)
-    );
-
-    hEqn.relax();
-    hEqn.solve(mesh.solver(h.select(finalIter)));
-
-    thermo.correct();
-
-    rad.correct();
-
-    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
-        << max(thermo.T()).value() << endl;
-
-}

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -1,6 +1,8 @@
     fvMesh& mesh = fluidRegions[i];
 
     rhoThermo& thermo = thermoFluid[i];
+    thermo.validate(args.executable(), "h", "e");
+
     volScalarField& rho = rhoFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
@@ -11,7 +13,6 @@
 
     volScalarField& p = thermo.p();
     const volScalarField& psi = thermo.psi();
-    volScalarField& h = thermo.he();
 
     volScalarField& p_rgh = p_rghFluid[i];
     const volScalarField& gh = ghFluid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H

@@ -9,8 +9,7 @@ if (oCorr == 0)
 }
 
 #include "UEqn.H"
-
+#include "EEqn.H"
-#include "hEqn.H"
 
 // --- PISO loop
 for (int corr=0; corr<nCorr; corr++)

applications/solvers/lagrangian/LTSReactingParcelFoam/UEqn.H

@@ -1,19 +0,0 @@
-    fvVectorMatrix UEqn
-    (
-        //pZones.ddt(rho, U)
-        fvm::ddt(rho, U)
-      + fvm::div(phi, U)
-      + turbulence->divDevRhoReff(U)
-     ==
-        rho.dimensionedInternalField()*g
-      + parcels.SU(U)
-    );
-
-    sources.constrain(UEqn);
-
-    pZones.addResistance(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve(UEqn == -fvc::grad(p) + sources(rho, U));
-    }

applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H

@@ -1,52 +0,0 @@
-tmp<fv::convectionScheme<scalar> > mvConvection
-(
-    fv::convectionScheme<scalar>::New
-    (
-        mesh,
-        fields,
-        phi,
-        mesh.divScheme("div(phi,Yi_h)")
-    )
-);
-
-combustion->correct();
-dQ = combustion->dQ();
-
-if (solveSpecies)
-{
-    label inertIndex = -1;
-    volScalarField Yt(0.0*Y[0]);
-
-    forAll(Y, i)
-    {
-        if (Y[i].name() != inertSpecie)
-        {
-            volScalarField& Yi = Y[i];
-
-            fvScalarMatrix YEqn
-            (
-                fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
-             ==
-                parcels.SYi(i, Yi)
-              + combustion->R(Yi)
-              + sources(rho, Yi)
-            );
-
-            sources.constrain(YEqn);
-
-            YEqn.solve(mesh.solver("Yi"));
-
-            Yi.max(0.0);
-            Yt += Yi;
-        }
-        else
-        {
-            inertIndex = i;
-        }
-    }
-
-    Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
-}

applications/solvers/lagrangian/LTSReactingParcelFoam/chemistry.H

@@ -1,28 +0,0 @@
-if (chemistry.chemistry())
-{
-    Info<< "Solving chemistry" << endl;
-
-    // update reaction rates
-    chemistry.calculate();
-
-    // turbulent time scale
-    if (turbulentReaction)
-    {
-        typedef DimensionedField<scalar, volMesh> dsfType;
-
-        const dimensionedScalar e0("e0", sqr(dimLength)/pow3(dimTime), SMALL);
-
-        const dsfType tk =
-            Cmix*sqrt(turbulence->muEff()/rho/(turbulence->epsilon() + e0));
-
-        const dsfType tc = chemistry.tc()().dimensionedInternalField();
-
-        kappa = tc/(tc + tk);
-    }
-    else
-    {
-        kappa = 1.0;
-    }
-
-    chemistrySh = kappa*chemistry.Sh()();
-}

applications/solvers/lagrangian/LTSReactingParcelFoam/createClouds.H

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

applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H

@@ -1,121 +0,0 @@
-    Info<< "Creating combustion model\n" << endl;
-
-    autoPtr<combustionModels::rhoCombustionModel> combustion
-    (
-        combustionModels::rhoCombustionModel::New(mesh)
-    );
-
-    rhoReactionThermo& thermo = combustion->thermo();
-
-    SLGThermo slgThermo(mesh, thermo);
-
-    basicMultiComponentMixture& composition = thermo.composition();
-    PtrList<volScalarField>& Y = composition.Y();
-
-    const word inertSpecie(thermo.lookup("inertSpecie"));
-
-    if (!composition.contains(inertSpecie))
-    {
-        FatalErrorIn(args.executable())
-            << "Specified inert specie '" << inertSpecie << "' not found in "
-            << "species list. Available species:" << composition.species()
-            << exit(FatalError);
-    }
-
-    volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.he();
-    const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.rho()
-    );
-
-    Info<< "\nReading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "compressibleCreatePhi.H"
-
-    dimensionedScalar rhoMax
-    (
-        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax")
-    );
-
-    dimensionedScalar rhoMin
-    (
-        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin")
-    );
-
-    Info<< "Creating turbulence model\n" << endl;
-    autoPtr<compressible::turbulenceModel> turbulence
-    (
-        compressible::turbulenceModel::New
-        (
-            rho,
-            U,
-            phi,
-            thermo
-        )
-    );
-
-    // Set the turbulence into the combustion model
-    combustion->setTurbulence(turbulence());
-
-    Info<< "Creating multi-variate interpolation scheme\n" << endl;
-    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-    forAll(Y, i)
-    {
-        fields.add(Y[i]);
-    }
-    fields.add(hs);
-
-    volScalarField dQ
-    (
-        IOobject
-        (
-            "dQ",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
-    );
-
-
-    volScalarField rDeltaT
-    (
-        IOobject
-        (
-            "rDeltaT",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("one", dimless/dimTime, 1),
-        zeroGradientFvPatchScalarField::typeName
-    );

applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H

@@ -1,25 +0,0 @@
-{
-    fvScalarMatrix hsEqn
-    (
-        fvm::ddt(rho, hs)
-      + mvConvection->fvmDiv(phi, hs)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
-     ==
-      - fvc::div(phi, 0.5*magSqr(U), "div(phi,K)")
-      + parcels.Sh(hs)
-      + radiation->Sh(thermo)
-      + combustion->Sh()
-      + sources(rho, hs)
-    );
-
-    sources.constrain(hsEqn);
-
-    hsEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-
-    Info<< "T gas min/max   = " << min(T).value() << ", "
-        << max(T).value() << endl;
-}

applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H

@@ -1,66 +0,0 @@
-{
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*(UEqn == sources(rho, U))().H();
-
-    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
-    if (pZones.size() == 0)
-    {
-        // ddtPhiCorr only used without porosity
-        phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
-    }
-
-    phiHbyA *= fvc::interpolate(rho);
-
-
-    fvScalarMatrix pDDtEqn
-    (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA)
-     ==
-        parcels.Srho()
-      + sources(psi, p, rho.name())
-    );
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            pDDtEqn
-          - fvm::laplacian(rho*rAU, p)
-        );
-
-        sources.constrain(pDDtEqn, rho.name());
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
-    }
-
-    p.relax();
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    sources.correct(U);
-
-    rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
-
-    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
-}

applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H

@@ -1,9 +0,0 @@
-dictionary additional = mesh.solutionDict().subDict("additional");
-
-// pressure work term for enthalpy equation
-bool pressureWork = additional.lookupOrDefault("pressureWork", true);
-bool pressureWorkTimeDerivative =
-    additional.lookupOrDefault("pressureWorkTimeDerivative", true);
-
-// flag to activate solve transport for each specie (Y vector)
-bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);

applications/solvers/lagrangian/LTSReactingParcelFoam/rhoEqn.H

@@ -1,50 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 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
-    rhoEqn
-
-Description
-    Solve the continuity for density.
-
-\*---------------------------------------------------------------------------*/
-
-{
-    fvScalarMatrix rhoEqn
-    (
-        fvm::ddt(rho)
-      + fvc::div(phi)
-      ==
-        parcels.Srho(rho)
-      + sources(rho)
-    );
-
-    sources.constrain(rhoEqn);
-
-    rhoEqn.solve();
-
-    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
-        << endl;
-}
-
-// ************************************************************************* //

applications/solvers/lagrangian/coalChemistryFoam/EEqn.H

@@ -0,0 +1,31 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        combustion->Sh()
+      + coalParcels.Sh(he)
+      + limestoneParcels.Sh(he)
+      + radiation->Sh(thermo)
+      + sources(rho, he)
+    );
+
+    EEqn.relax();
+    sources.constrain(EEqn);
+    EEqn.solve();
+
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hsEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -6,6 +6,7 @@
     );
 
     psiReactionThermo& thermo = combustion->thermo();
+    thermo.validate(args.executable(), "h", "e");
 
     SLGThermo slgThermo(mesh, thermo);
 
@@ -23,7 +24,6 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.he();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -33,7 +33,7 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(hs);
+    fields.add(thermo.he());
 
     volScalarField rho
     (

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -1,29 +0,0 @@
-{
-    fvScalarMatrix hsEqn
-    (
-        fvm::ddt(rho, hs)
-      + mvConvection->fvmDiv(phi, hs)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + combustion->Sh()
-      + coalParcels.Sh(hs)
-      + limestoneParcels.Sh(hs)
-      + radiation->Sh(thermo)
-      + sources(rho, hs)
-    );
-
-    hsEqn.relax();
-
-    sources.constrain(hsEqn);
-
-    hsEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-
-    Info<< "T gas min/max   = " << min(T).value() << ", "
-        << max(T).value() << endl;
-}

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Allwmake

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+wmake icoUncoupledKinematicParcelDyMFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I../icoUncoupledKinematicParcelFoam \
+    -I.. \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/files

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

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/Make/options

@@ -1,54 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I${LIB_SRC}/meshTools/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
-    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(LIB_SRC)/fieldSources/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam
-
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lcompressibleTurbulenceModel \
-    -lcompressibleRASModels \
-    -lcompressibleLESModels \
-    -llagrangian \
-    -llagrangianIntermediate \
-    -lspecie \
-    -lfluidThermophysicalModels \
-    -lliquidProperties \
-    -lliquidMixtureProperties \
-    -lsolidProperties \
-    -lsolidMixtureProperties \
-    -lthermophysicalFunctions \
-    -lreactionThermophysicalModels \
-    -lSLGThermo \
-    -lchemistryModel \
-    -lradiationModels \
-    -lODE \
-    -lregionModels \
-    -lsurfaceFilmModels \
-    -lcombustionModels \
-    -lfieldSources \
-    -lsampling
-

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -1,55 +0,0 @@
-
-tmp<fv::convectionScheme<scalar> > mvConvection
-(
-    fv::convectionScheme<scalar>::New
-    (
-        mesh,
-        fields,
-        phi,
-        mesh.divScheme("div(phi,Yi_h)")
-    )
-);
-
-combustion->correct();
-dQ = combustion->dQ();
-
-if (solveSpecies)
-{
-    label inertIndex = -1;
-    volScalarField Yt(0.0*Y[0]);
-
-    forAll(Y, i)
-    {
-        if (Y[i].name() != inertSpecie)
-        {
-            volScalarField& Yi = Y[i];
-
-            fvScalarMatrix YiEqn
-            (
-                fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(turbulence->muEff(), Yi)
-              ==
-                parcels.SYi(i, Yi)
-              + combustion->R(Yi)
-              + sources(rho, Yi)
-            );
-
-            YiEqn.relax();
-
-            sources.constrain(YiEqn);
-
-            YiEqn.solve(mesh.solver("Yi"));
-
-            Yi.max(0.0);
-            Yt += Yi;
-        }
-        else
-        {
-            inertIndex = i;
-        }
-    }
-
-    Y[inertIndex] = scalar(1) - Yt;
-    Y[inertIndex].max(0.0);
-}

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createClouds.H

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

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createExplicitSources.H

@@ -1,2 +0,0 @@
-Info<< "Creating sources\n" << endl;
-IObasicSourceList sources(mesh);

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createPorousZones.H

@@ -1,3 +0,0 @@
-    Info<< "Creating porous zones" << nl << endl;
-
-    porousZones pZones(mesh);

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H

@@ -1,32 +0,0 @@
-{
-    fvScalarMatrix hsEqn
-    (
-        fvm::ddt(rho, hs)
-      + mvConvection->fvmDiv(phi, hs)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
-     ==
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + parcels.Sh(hs)
-      + radiation->Sh(thermo)
-      + combustion->Sh()
-      + sources(rho, hs)
-    );
-
-    if (pressureWorkTimeDerivative)
-    {
-        hsEqn -= dpdt;
-    }
-
-    hsEqn.relax();
-
-    sources.constrain(hsEqn);
-
-    hsEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-
-    Info<< "T gas min/max   = " << min(T).value() << ", "
-        << max(T).value() << endl;
-}

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -1,66 +0,0 @@
-{
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    volVectorField HbyA("HbyA", U);
-    HbyA = rAU*(UEqn == sources(rho, U))().H();
-
-    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
-    if (pZones.size() == 0)
-    {
-        // ddtPhiCorr only used without porosity
-        phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
-    }
-
-    phiHbyA *= fvc::interpolate(rho);
-
-
-    fvScalarMatrix pDDtEqn
-    (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA)
-     ==
-        parcels.Srho()
-      + sources(psi, p, rho.name())
-    );
-
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix pEqn
-        (
-            pDDtEqn
-          - fvm::laplacian(rho*rAU, p)
-        );
-
-        sources.constrain(pDDtEqn, rho.name());
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
-    }
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    sources.correct(U);
-
-    K = 0.5*magSqr(U);
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-}

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C

@@ -1,126 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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
-    porousExplicitSourceReactingParcelFoam
-
-Description
-    Transient PIMPLE solver for compressible, laminar or turbulent flow with
-    reacting multiphase Lagrangian parcels for porous media, including explicit
-    sources for mass, momentum and energy
-
-    The solver includes:
-    - reacting multiphase parcel cloud
-    - porous media
-    - mass, momentum and energy sources
-
-    Note: ddtPhiCorr not used here when porous zones are active
-    - not well defined for porous calculations
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "turbulenceModel.H"
-#include "basicReactingMultiphaseCloud.H"
-#include "rhoCombustionModel.H"
-#include "radiationModel.H"
-#include "porousZones.H"
-#include "IObasicSourceList.H"
-#include "SLGThermo.H"
-#include "pimpleControl.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
-    #include "createFields.H"
-    #include "createRadiationModel.H"
-    #include "createClouds.H"
-    #include "createExplicitSources.H"
-    #include "createPorousZones.H"
-    #include "initContinuityErrs.H"
-    #include "readTimeControls.H"
-    #include "compressibleCourantNo.H"
-    #include "setInitialDeltaT.H"
-
-    pimpleControl pimple(mesh);
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        #include "readTimeControls.H"
-        #include "readAdditionalSolutionControls.H"
-        #include "compressibleCourantNo.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        parcels.evolve();
-
-        #include "rhoEqn.H"
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            #include "UEqn.H"
-            #include "YEqn.H"
-            #include "hsEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.turbCorr())
-            {
-                turbulence->correct();
-            }
-        }
-
-        rho = thermo.rho();
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H

@@ -0,0 +1,29 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        parcels.Sh(he)
+      + surfaceFilm.Sh()
+      + radiation->Sh(thermo)
+      + combustion->Sh()
+    );
+
+    EEqn.relax();
+    EEqn.solve();
+
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -6,6 +6,7 @@
     );
 
     psiReactionThermo& thermo = combustion->thermo();
+    thermo.validate(args.executable(), "h", "e");
 
     SLGThermo slgThermo(mesh, thermo);
 
@@ -29,7 +30,6 @@
     );
 
     volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.he();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -107,7 +107,7 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(hs);
+    fields.add(thermo.he());
 
     IOdictionary additionalControlsDict
     (

applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H

@@ -1,25 +0,0 @@
-{
-    fvScalarMatrix hsEqn
-    (
-        fvm::ddt(rho, hs)
-      + mvConvection->fvmDiv(phi, hs)
-      - fvm::laplacian(turbulence->alphaEff(), hs)
-     ==
-        dpdt
-      - (fvc::ddt(rho, K) + fvc::div(phi, K))
-      + parcels.Sh(hs)
-      + surfaceFilm.Sh()
-      + radiation->Sh(thermo)
-      + combustion->Sh()
-    );
-
-    hsEqn.relax();
-
-    hsEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-
-    Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
-}

applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C

@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
             {
                 #include "UEqn.H"
                 #include "YEqn.H"
-                #include "hsEqn.H"
+                #include "EEqn.H"
 
                 // --- Pressure corrector loop
                 while (pimple.correct())

applications/solvers/lagrangian/reactingParcelFoam/Allwmake

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+wmake LTSReactingParcelFoam
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/lagrangian/reactingParcelFoam/EEqn.H

@@ -0,0 +1,30 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ep", p/rho))
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        parcels.Sh(he)
+      + radiation->Sh(thermo)
+      + combustion->Sh()
+      + sources(rho, he)
+    );
+
+    EEqn.relax();
+    sources.constrain(EEqn);
+    EEqn.solve();
+
+    thermo.correct();
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
     #include "readTimeControls.H"
     #include "readAdditionalSolutionControls.H"
     #include "createFields.H"
+    #include "createRDeltaT.H"
     #include "createRadiationModel.H"
     #include "createClouds.H"
     #include "createExplicitSources.H"
@@ -93,7 +94,7 @@ int main(int argc, char *argv[])
 
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hsEqn.H"
+            #include "EEqn.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I.. \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/createRDeltaT.H

@@ -0,0 +1,14 @@
+    volScalarField rDeltaT
+    (
+        IOobject
+        (
+            "rDeltaT",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless/dimTime, 1),
+        zeroGradientFvPatchScalarField::typeName
+    );

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/readTimeControls.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/timeScales.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/lagrangian/reactingParcelFoam/Make/options

@@ -1,10 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
-    -I${LIB_SRC}/sampling/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
@@ -22,13 +22,13 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
-    -lsampling \
     -lcompressibleTurbulenceModel \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
@@ -48,5 +48,6 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
+    -lcombustionModels \
     -lfieldSources \
-    -lcombustionModels
+    -lsampling

applications/solvers/lagrangian/reactingParcelFoam/UEqn.H

@@ -1,5 +1,6 @@
     fvVectorMatrix UEqn
     (
+        //pZones.ddt(rho, U)
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
@@ -12,6 +13,8 @@
 
     sources.constrain(UEqn);
 
+    pZones.addResistance(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p) + sources(rho, U));

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -5,14 +5,15 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,Yi_hs)")
+        mesh.divScheme("div(phi,Yi_h)")
     )
 );
 
+combustion->correct();
+dQ = combustion->dQ();
 
+if (solveSpecies)
 {
-    combustion->correct();
-    dQ = combustion->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -22,7 +23,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         {
             volScalarField& Yi = Y[i];
 
-            fvScalarMatrix YiEqn
+            fvScalarMatrix YEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
@@ -33,11 +34,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               + sources(rho, Yi)
             );
 
-            YiEqn.relax();
+            YEqn.relax();
 
-            sources.constrain(YiEqn);
+            sources.constrain(YEqn);
 
-            YiEqn.solve(mesh.solver("Yi"));
+            YEqn.solve(mesh.solver("Yi"));
 
             Yi.max(0.0);
             Yt += Yi;

applications/solvers/lagrangian/reactingParcelFoam/createClouds.H

@@ -1,5 +1,5 @@
 Info<< "\nConstructing reacting cloud" << endl;
-basicReactingCloud parcels
+basicReactingMultiphaseCloud parcels
 (
     "reactingCloud1",
     rho,

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -1,11 +1,12 @@
     Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<combustionModels::psiCombustionModel> combustion
+    autoPtr<combustionModels::rhoCombustionModel> combustion
     (
-        combustionModels::psiCombustionModel::New(mesh)
+        combustionModels::rhoCombustionModel::New(mesh)
     );
 
-    psiReactionThermo& thermo = combustion->thermo();
+    rhoReactionThermo& thermo = combustion->thermo();
+    thermo.validate(args.executable(), "h", "e");
 
     SLGThermo slgThermo(mesh, thermo);
 
@@ -23,7 +24,6 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& hs = thermo.he();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -56,6 +56,22 @@
 
     #include "compressibleCreatePhi.H"
 
+    dimensionedScalar rhoMax
+    (
+        "rhoMax",
+        dimDensity,
+        mesh.solutionDict().subDict("PIMPLE")
+       .lookupOrDefault<scalar>("rhoMax", GREAT)
+    );
+
+    dimensionedScalar rhoMin
+    (
+        "rhoMin",
+        dimDensity,
+        mesh.solutionDict().subDict("PIMPLE")
+       .lookupOrDefault<scalar>("rhoMin", -GREAT)
+    );
+
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -86,13 +102,15 @@
 
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "Creating multi-variate interpolation scheme\n" << endl;
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll(Y, i)
     {
         fields.add(Y[i]);
     }
-    fields.add(hs);
+    fields.add(thermo.he());
 
     volScalarField dQ
     (
@@ -107,6 +125,3 @@
         mesh,
         dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
-
-    Info<< "Creating sources\n" << endl;
-    IObasicSourceList sources(mesh);