diff --git a/EEqn.H b/EEqn.H
index 9267c9a..236a797 100644
--- a/EEqn.H
+++ b/EEqn.H
@@ -1,23 +1,62 @@
 {
+    typedef multiComponentMixture<gasHThermoPhysics> MMix;
+
+    MMix &janafComposition = dynamic_cast<MMix&>(composition);
+
     volScalarField& he = thermo.he();
 
+    volScalarField hsi(he);
+
+    volVectorField hDiffusionSrc("rho*Sum(hYV)", he * Vc * 0.0);
+
+    forAll(Y, i)
+    {
+        const MMix::thermoType &tti = janafComposition.getLocalThermo(i);
+
+        forAll(hsi, cellI)
+        {
+            const scalar pI = p[cellI];
+            const scalar TI = T[cellI];
+
+            hsi[cellI] = tti.HE(pI,TI);
+        }
+
+        forAll(hsi.boundaryFieldRef(), patchI)
+        {
+            volScalarField::Patch &hsiP = hsi.boundaryFieldRef()[patchI];
+            const volScalarField::Patch &pP = p.boundaryField()[patchI];
+            const volScalarField::Patch &TP = T.boundaryField()[patchI];
+
+            forAll(hsiP, faceI)
+            {
+                const scalar pI = pP[faceI];
+                const scalar TI = TP[faceI];
+
+                hsiP[faceI] = tti.HE(pI,TI);
+            }
+        }
+
+        hDiffusionSrc += hsi * (- diff.D(i) * fvc::grad(Y[i]) + Y[i] * Vc);
+    }
+
+    tmp<volVectorField> tTauU(diff.mu()*(U&(Foam::dev2(Foam::T(gradU)) + gradU)));
+
+    hDiffusionSrc *= rho;
+
+    volVectorField hWorkSrc1(rho * U);
+
+    volVectorField hWorkSrc2(rho * Vc);
+
     fvScalarMatrix EEqn
     (
         fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
       + fvc::ddt(rho, K) + fvc::div(phi, K)
-      + (
-            he.name() == "e"
-          ? fvc::div
-            (
-                fvc::absolute(phi/fvc::interpolate(rho), U),
-                p,
-                "div(phiv,p)"
-            )
-          : -dpdt
-        )
-      - fvm::laplacian(turbulence->alphaEff(), he)
+      - dpdt
+      - fvc::laplacian(diff.k(), T)
+      + fvc::div(hDiffusionSrc)
      ==
         reaction->Sh()
+      + fvc::div(tTauU)
       + fvOptions(rho, he)
     );
 
diff --git a/UEqn.H b/UEqn.H
index 849fe89..1058faa 100644
--- a/UEqn.H
+++ b/UEqn.H
@@ -2,13 +2,14 @@
 
 MRF.correctBoundaryVelocity(U);
 
+volTensorField gradU(fvc::grad(U));
 
 tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
   + (
-     - fvc::div(diff.mu()*Foam::dev2(Foam::T(fvc::grad(U))))
+     - fvc::div(diff.mu()*Foam::dev2(Foam::T(gradU)))
      - fvm::laplacian(diff.mu(), U)
     )
  ==
diff --git a/createFields.H b/createFields.H
index a2cb408..d6f13dc 100644
--- a/createFields.H
+++ b/createFields.H
@@ -6,7 +6,7 @@ autoPtr<combustionModels::psiCombustionModel> reaction
 );
 
 psiReactionThermo& thermo = reaction->thermo();
-thermo.validate(args.executable(), "h", "e");
+thermo.validate(args.executable(), "h");
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index 3c9176d..39c9e98 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -298,7 +298,7 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
             IOobject::AUTO_WRITE
         ),
         thermo_.composition().Y(0).mesh(),
-        dimensionedScalar("zero", dimPressure*dimTime, 0.0)
+        dimensionedScalar("zero", dimDynamicViscosity, 0.0)
     ),
     k_
     (
@@ -311,7 +311,7 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
             IOobject::AUTO_WRITE
         ),
         thermo_.composition().Y(0).mesh(),
-        dimensionedScalar("zero", dimArea/dimTime, 0.0)
+        dimensionedScalar("zero", dimForce/dimTime/dimTemperature, 0.0)
     ),
     neutrals_(thermo_.composition().species().size()),
     ions_(thermo_.composition().species().size()),
@@ -673,7 +673,7 @@ void Foam::diffusivityModel::correct()
         }
 
         mu_[celli] = mixAvgMu(muI, localX, Wpure);
-        k_[celli] = mixAvgK(kI, localX) / rhoi / Cpi;
+        k_[celli] = mixAvgK(kI, localX);
     }
 
 
@@ -750,7 +750,7 @@ void Foam::diffusivityModel::correct()
             }
 
             mu_.boundaryFieldRef()[patchi][facei] = mixAvgMu(muI, localX, Wpure);
-            k_.boundaryFieldRef()[patchi][facei] = mixAvgK(kI, localX) / rhoi / Cpi;
+            k_.boundaryFieldRef()[patchi][facei] = mixAvgK(kI, localX);
         }
     }
 
diff --git a/eReactingFoam.C b/eReactingFoam.C
index 584372d..7234f2a 100644
--- a/eReactingFoam.C
+++ b/eReactingFoam.C
@@ -38,6 +38,9 @@ Description
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
+#include "thermoPhysicsTypes.H"
+#include "multiComponentMixture.H"
+
 #include "diffusivityModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //