Thermodynamics: Updated combustion solvers to use EEqn and support h/e or ha/ea

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,7 +5,7 @@
         psiuReactionThermo::New(mesh)
     );
     psiuReactionThermo& thermo = pThermo();
-    thermo.validate(args.executable(), "ha");
+    thermo.validate(args.executable(), "ha", "ea");
 
     basicMultiComponentMixture& composition = thermo.composition();
 
@@ -24,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
     (
@@ -198,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,7 +5,7 @@
         psiuReactionThermo::New(mesh)
     );
     psiuReactionThermo& thermo = pThermo();
-    thermo.validate(args.executable(), "ha");
+    thermo.validate(args.executable(), "ha", "ea");
 
     basicMultiComponentMixture& composition = thermo.composition();
 
@@ -24,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
@@ -140,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/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

@@ -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())

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/Allclean

@@ -6,6 +6,7 @@ cd ${0%/*} || exit 1    # run from this directory
 cleanCase
 
 rm -rf VTK
+rm -rf 0
 rm -rf constant/polyMesh/sets constant/polyMesh/faces
 rm -rf constant/polyMesh/faces
 rm -rf constant/polyMesh/points

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/Allrun

@@ -5,11 +5,13 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
+rm -rf 0
+cp -r 0.org 0
 runApplication blockMesh
 runApplication changeDictionary
 runApplication topoSet
 
-runApplication PDRMesh
+runApplication PDRMesh -overwrite
 
 # Run
 runApplication `getApplication`

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/polyMesh/boundary

@@ -44,6 +44,7 @@ FoamFile
     baffleCyclic_half0
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          0;
         startFace       24810;
         matchTolerance  0.0001;
@@ -52,6 +53,7 @@ FoamFile
     baffleCyclic_half1
     {
         type            cyclic;
+        inGroups        1(cyclic);
         nFaces          0;
         startFace       24810;
         matchTolerance  0.0001;

tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/polyMesh/boundary

@@ -38,6 +38,7 @@ FoamFile
     frontAndBack
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          45000;
         startFace       45284;
     }