Corrected handling of density.

applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H

@@ -9,7 +9,7 @@
       - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
-    pZones.addEnthalpySource(thermo, hEqn);
+    pZones.addEnthalpySource(thermo, rho, hEqn);
 
     hEqn.relax();
 

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C

@@ -73,6 +73,7 @@ Foam::thermalPorousZone::thermalPorousZone
 void Foam::thermalPorousZone::addEnthalpySource
 (
     const basicThermo& thermo,
+    const volScalarField& rho,
     fvScalarMatrix& hEqn
 ) const
 {
@@ -85,7 +86,6 @@ void Foam::thermalPorousZone::addEnthalpySource
     const scalarField& V = mesh().V();
     scalarField& hDiag = hEqn.diag();
     scalarField& hSource = hEqn.source();
-    const scalarField& rho = thermo.rho();
 
     scalarField hZone = thermo.h(scalarField(cells.size(), T_.value()), cells);
     scalar rate = 1e6;

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H

@@ -138,6 +138,7 @@ public:
         void addEnthalpySource
         (
             const basicThermo& thermo,
+            const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;
 };

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C

@@ -50,12 +50,13 @@ Foam::thermalPorousZones::thermalPorousZones
 void Foam::thermalPorousZones::addEnthalpySource
 (
     const basicThermo& thermo,
+    const volScalarField& rho,
     fvScalarMatrix& hEqn
 ) const
 {
     forAll(*this, i)
     {
-        operator[](i).addEnthalpySource(thermo, hEqn);
+        operator[](i).addEnthalpySource(thermo, rho, hEqn);
     }
 }
 

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H

@@ -92,6 +92,7 @@ public:
         void addEnthalpySource
         (
             const basicThermo& thermo,
+            const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;
 };