From 40f141e2c96af2d60355e76daaa83747d77014d2 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Fri, 12 Feb 2010 17:05:15 +0000
Subject: [PATCH] ENH: Updated enthalpy equations for solvers with chemistry
 using updated      sensible enthalpy thermo packages

- Enthalpy source term now retrieved from the chemistry model (scaled by
  kappa for the PaSR model)
---
 .../dieselEngineFoam/createFields.H           | 20 +++++-
 .../dieselEngineFoam/dieselEngineFoam.C       |  6 +-
 .../combustion/dieselEngineFoam/hEqn.H        | 13 ----
 .../combustion/dieselEngineFoam/hsEqn.H       | 14 ++++
 .../combustion/dieselFoam/dieselFoam.C        |  2 +-
 .../combustion/reactingFoam/Make/options      |  1 -
 .../solvers/combustion/reactingFoam/UEqn.H    | 15 ++++
 .../combustion/reactingFoam/chemistry.H       |  2 +
 .../combustion/reactingFoam/createFields.H    | 21 ++++--
 .../solvers/combustion/reactingFoam/hsEqn.H   | 20 ++++++
 .../solvers/combustion/reactingFoam/pEqn.H    | 72 +++++++++++++++++++
 .../combustion/reactingFoam/reactingFoam.C    |  4 +-
 .../solvers/combustion/reactingFoam/rhoEqn.H  | 42 +++++++++++
 .../combustion/rhoReactingFoam/Make/options   |  1 -
 .../solvers/combustion/rhoReactingFoam/UEqn.H | 15 ++++
 .../combustion/rhoReactingFoam/chemistry.H    |  2 +
 .../combustion/rhoReactingFoam/createFields.H | 20 +++++-
 .../combustion/rhoReactingFoam/hsEqn.H        | 19 +++++
 .../rhoReactingFoam/rhoReactingFoam.C         |  2 +-
 .../lagrangian/coalChemistryFoam/chemistry.H  |  2 +
 .../coalChemistryFoam/coalChemistryFoam.C     | 13 +---
 .../coalChemistryFoam/createFields.H          | 22 ++++--
 .../lagrangian/coalChemistryFoam/hEqn.H       | 22 ------
 .../lagrangian/coalChemistryFoam/hsEqn.H      | 26 +++++++
 .../chemistry.H                               |  2 +
 .../createFields.H                            | 20 +++++-
 .../{hEqn.H => hsEqn.H}                       |  9 +--
 .../porousExplicitSourceReactingParcelFoam.C  |  6 +-
 .../lagrangian/reactingParcelFoam/chemistry.H |  2 +
 .../reactingParcelFoam/createFields.H         | 20 +++++-
 .../lagrangian/reactingParcelFoam/hEqn.H      | 20 ------
 .../lagrangian/reactingParcelFoam/hsEqn.H     | 24 +++++++
 .../reactingParcelFoam/reactingParcelFoam.C   |  7 +-
 33 files changed, 376 insertions(+), 110 deletions(-)
 delete mode 100644 applications/solvers/combustion/dieselEngineFoam/hEqn.H
 create mode 100644 applications/solvers/combustion/dieselEngineFoam/hsEqn.H
 create mode 100644 applications/solvers/combustion/reactingFoam/UEqn.H
 create mode 100644 applications/solvers/combustion/reactingFoam/hsEqn.H
 create mode 100644 applications/solvers/combustion/reactingFoam/pEqn.H
 create mode 100644 applications/solvers/combustion/reactingFoam/rhoEqn.H
 create mode 100644 applications/solvers/combustion/rhoReactingFoam/UEqn.H
 create mode 100644 applications/solvers/combustion/rhoReactingFoam/hsEqn.H
 delete mode 100644 applications/solvers/lagrangian/coalChemistryFoam/hEqn.H
 create mode 100644 applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
 rename applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/{hEqn.H => hsEqn.H} (83%)
 delete mode 100644 applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
 create mode 100644 applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H

diff --git a/applications/solvers/combustion/dieselEngineFoam/createFields.H b/applications/solvers/combustion/dieselEngineFoam/createFields.H
index a75c6a7d07a..37aeca4c091 100644
--- a/applications/solvers/combustion/dieselEngineFoam/createFields.H
+++ b/applications/solvers/combustion/dieselEngineFoam/createFields.H
@@ -6,7 +6,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -50,7 +50,7 @@ volVectorField U
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
 
 
 #include "compressibleCreatePhi.H"
@@ -92,4 +92,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
+
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);
diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
index 6507c117a01..2e6a9c4b05c 100644
--- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
+++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    dieselFoam
+    dieselEngineFoam
 
 Description
     Solver for diesel engine spray and combustion.
@@ -103,13 +103,15 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
         }
 
+        chemistrySh = kappa*chemistry.Sh()();
+
         #include "rhoEqn.H"
         #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/combustion/dieselEngineFoam/hEqn.H b/applications/solvers/combustion/dieselEngineFoam/hEqn.H
deleted file mode 100644
index 0406f7fbd49..00000000000
--- a/applications/solvers/combustion/dieselEngineFoam/hEqn.H
+++ /dev/null
@@ -1,13 +0,0 @@
-{
-    solve
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-       DpDt
-     + dieselSpray.heatTransferSource()
-    );
-
-    thermo.correct();
-}
diff --git a/applications/solvers/combustion/dieselEngineFoam/hsEqn.H b/applications/solvers/combustion/dieselEngineFoam/hsEqn.H
new file mode 100644
index 00000000000..7ae59feb819
--- /dev/null
+++ b/applications/solvers/combustion/dieselEngineFoam/hsEqn.H
@@ -0,0 +1,14 @@
+{
+    solve
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+       DpDt
+     + dieselSpray.heatTransferSource()().dimensionedInternalField()
+     + chemistrySh
+    );
+
+    thermo.correct();
+}
diff --git a/applications/solvers/combustion/dieselFoam/dieselFoam.C b/applications/solvers/combustion/dieselFoam/dieselFoam.C
index 4769a00c77e..1e6ddfe0c19 100644
--- a/applications/solvers/combustion/dieselFoam/dieselFoam.C
+++ b/applications/solvers/combustion/dieselFoam/dieselFoam.C
@@ -100,7 +100,7 @@ int main(int argc, char *argv[])
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/combustion/reactingFoam/Make/options b/applications/solvers/combustion/reactingFoam/Make/options
index e982f635a14..d85350297af 100644
--- a/applications/solvers/combustion/reactingFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
diff --git a/applications/solvers/combustion/reactingFoam/UEqn.H b/applications/solvers/combustion/reactingFoam/UEqn.H
new file mode 100644
index 00000000000..9697c6e1ed2
--- /dev/null
+++ b/applications/solvers/combustion/reactingFoam/UEqn.H
@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }
diff --git a/applications/solvers/combustion/reactingFoam/chemistry.H b/applications/solvers/combustion/reactingFoam/chemistry.H
index a1a978210dd..8bfa872d7b3 100644
--- a/applications/solvers/combustion/reactingFoam/chemistry.H
+++ b/applications/solvers/combustion/reactingFoam/chemistry.H
@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 0b2238031ca..bd8f7ea1fb1 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -5,7 +5,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,8 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
-
+volScalarField& hs = thermo.hs();
+const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
@@ -81,5 +81,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);
diff --git a/applications/solvers/combustion/reactingFoam/hsEqn.H b/applications/solvers/combustion/reactingFoam/hsEqn.H
new file mode 100644
index 00000000000..d1904badee5
--- /dev/null
+++ b/applications/solvers/combustion/reactingFoam/hsEqn.H
@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+//      - fvm::laplacian(turbulence->alphaEff(), hs)
+      - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
+     ==
+        DpDt
+      + chemistrySh
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
new file mode 100644
index 00000000000..62dbd45b12d
--- /dev/null
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -0,0 +1,72 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+//          ==
+//            RRTot
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+//          ==
+//            RRTot
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index 2d4ae7589ad..e9fe757e3ba 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.C
@@ -73,9 +73,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-
-            #define Db turbulence->alphaEff()
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/combustion/reactingFoam/rhoEqn.H b/applications/solvers/combustion/reactingFoam/rhoEqn.H
new file mode 100644
index 00000000000..e7ef140dc9a
--- /dev/null
+++ b/applications/solvers/combustion/reactingFoam/rhoEqn.H
@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    solve
+    (
+        fvm::ddt(rho) + fvc::div(phi)
+//      ==
+//        RRTot
+    );
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/combustion/rhoReactingFoam/Make/options b/applications/solvers/combustion/rhoReactingFoam/Make/options
index e982f635a14..d85350297af 100644
--- a/applications/solvers/combustion/rhoReactingFoam/Make/options
+++ b/applications/solvers/combustion/rhoReactingFoam/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
diff --git a/applications/solvers/combustion/rhoReactingFoam/UEqn.H b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
new file mode 100644
index 00000000000..9697c6e1ed2
--- /dev/null
+++ b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }
diff --git a/applications/solvers/combustion/rhoReactingFoam/chemistry.H b/applications/solvers/combustion/rhoReactingFoam/chemistry.H
index a1a978210dd..8bfa872d7b3 100644
--- a/applications/solvers/combustion/rhoReactingFoam/chemistry.H
+++ b/applications/solvers/combustion/rhoReactingFoam/chemistry.H
@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/combustion/rhoReactingFoam/createFields.H b/applications/solvers/combustion/rhoReactingFoam/createFields.H
index d2f6f5e4d49..d4119b8c926 100644
--- a/applications/solvers/combustion/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H
@@ -5,7 +5,7 @@ autoPtr<rhoChemistryModel> pChemistry
 );
 rhoChemistryModel& chemistry = pChemistry();
 
-hReactionThermo& thermo = chemistry.thermo();
+hsReactionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,7 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
+const volScalarField& T = thermo.T();
 
 
 #include "compressibleCreatePhi.H"
@@ -81,5 +82,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);
diff --git a/applications/solvers/combustion/rhoReactingFoam/hsEqn.H b/applications/solvers/combustion/rhoReactingFoam/hsEqn.H
new file mode 100644
index 00000000000..81bd6ea9d7b
--- /dev/null
+++ b/applications/solvers/combustion/rhoReactingFoam/hsEqn.H
@@ -0,0 +1,19 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+      + chemistrySh
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index cc37dd09c18..79b76e014d6 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H b/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
index 5489dfefbaa..50fc7f575b5 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/chemistry.H
@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 389a0659350..872664f3a52 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -84,16 +84,8 @@ int main(int argc, char *argv[])
 
         coalParcels.evolve();
 
-        coalParcels.info();
-
-        Info<< endl;
-
         limestoneParcels.evolve();
 
-        limestoneParcels.info();
-
-        Info<< endl;
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -102,16 +94,13 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
index d0a284a0e17..65be09775ca 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/createFields.H
@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -32,7 +32,7 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
 
     volScalarField rho
     (
@@ -133,5 +133,19 @@
         "energy",
         mesh,
         dimEnergy/dimTime/dimVolume,
-        "h"
+        "hs"
+    );
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
     );
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hEqn.H
deleted file mode 100644
index 2eaf6315cfb..00000000000
--- a/applications/solvers/lagrangian/coalChemistryFoam/hEqn.H
+++ /dev/null
@@ -1,22 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  coalParcels.Sh()
-     +  limestoneParcels.Sh()
-     +  enthalpySource.Su()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
new file mode 100644
index 00000000000..ced6b640044
--- /dev/null
+++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
@@ -0,0 +1,26 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  coalParcels.Sh()
+     +  limestoneParcels.Sh()
+     +  enthalpySource.Su()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hsEqn.relax();
+
+    hsEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
index 5489dfefbaa..50fc7f575b5 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H
@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
index 68c0f8aee8c..cf3f484bdbc 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H
@@ -6,7 +6,7 @@
     );
     rhoChemistryModel& chemistry = pChemistry();
 
-    hReactionThermo& thermo = chemistry.thermo();
+    hsReactionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -88,4 +88,18 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    );
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
similarity index 83%
rename from applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H
rename to applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
index 0715989df7a..73c276c2942 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
@@ -32,14 +32,15 @@
     {
         solve
         (
-            fvm::ddt(rho, h)
-          + mvConvection->fvmDiv(phi, h)
-          - fvm::laplacian(turbulence->alphaEff(), h)
+            fvm::ddt(rho, hs)
+          + mvConvection->fvmDiv(phi, hs)
+          - fvm::laplacian(turbulence->alphaEff(), hs)
          ==
             pWork()
           + parcels.Sh()
-          + radiation->Sh(thermo)
+          + radiation->Shs(thermo)
           + energySource.Su()
+          + chemistrySh
         );
 
         thermo.correct();
diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
index b345ffe28a0..436fda7a8bb 100644
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C
@@ -33,7 +33,7 @@ Description
     The solver includes:
     - reacting multiphase parcel cloud
     - porous media
-    - point mass sources
+    - mass, momentum and energy sources
     - polynomial based, incompressible thermodynamics (f(T))
 
     Note: ddtPhiCorr not used here when porous zones are active
@@ -89,13 +89,11 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
         #include "UEqn.H"
         #include "YEqn.H"
-        #include "hEqn.H"
+        #include "hsEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
index 5489dfefbaa..50fc7f575b5 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/chemistry.H
@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
index e825c2cffbc..e071ac8b3d3 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/createFields.H
@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -94,4 +94,18 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistry::Sh", dimEnergy/dimTime/dimVolume, 0.0)
+    );
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
deleted file mode 100644
index d686e452f46..00000000000
--- a/applications/solvers/lagrangian/reactingParcelFoam/hEqn.H
+++ /dev/null
@@ -1,20 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  parcels.Sh()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
new file mode 100644
index 00000000000..ce12ec3ad48
--- /dev/null
+++ b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
@@ -0,0 +1,24 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  parcels.Sh()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
index 4dfd367d092..ad3b5b461c4 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C
@@ -74,8 +74,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -88,12 +86,9 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-                #include "hEqn.H"
+                #include "hsEqn.H"
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();
-- 
GitLab