diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H
index bbbd079ef2..3ce6b03b9d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H
@@ -15,6 +15,9 @@
     const volScalarField& kappa = tkappa();
     //const volSymmTensorField& K = tK();
 
+    tmp<volScalarField> trhoCp = cp*rho;
+    const volScalarField& rhoCp = trhoCp();
+
     volScalarField& T = thermo.T();
 
     IObasicSourceList& sources = solidHeatSources[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H
index 2dd09b974f..5d4bb47b22 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H
@@ -4,7 +4,7 @@
         tmp<fvScalarMatrix> TEqn
         (
           - fvm::laplacian(betav*kappa, T, "laplacian(K,T)")
-          + sources(rho, T)
+          + sources(rhoCp, T)
         );
 
         TEqn().relax();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H
index 6f4834000d..2ade49d0ac 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H
@@ -15,6 +15,9 @@
     const volScalarField& kappa = tkappa();
     //const volSymmTensorField& K = tK();
 
+    tmp<volScalarField> trhoCp = cp*rho;
+    const volScalarField& rhoCp = trhoCp();
+
     volScalarField& T = thermo.T();
 
     IObasicSourceList& sources = solidHeatSources[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H
index e8941cf2ad..dfe0bff675 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H
@@ -10,7 +10,7 @@ if (finalIter)
         (
             fvm::ddt(betav*rho*cp, T)
           - fvm::laplacian(betav*kappa, T, "laplacian(K,T)")
-          + sources(rho, T)
+          + sources(rhoCp, T)
         );
 
         TEqn().relax();