updates for new thermo framework

applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C

@@ -27,13 +27,13 @@ Application
 
 Description
     Transient Solver for buoyant, turbulent flow of compressible fluids for
-    ventilation and heat-transfer.  Turbulence is modelled using a run-time
-    selectable compressible RAS model.
+    ventilation and heat-transfer. Turbulence is modelled using a run-time
+    selectable compressible RAS or LES model.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 
@@ -81,6 +81,8 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
+        rho = thermo->rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H

     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicRhoThermo> thermo
     (
-        basicThermo::New(mesh)
+        basicRhoThermo::New(mesh)
     );
 
     volScalarField rho
@@ -62,3 +62,6 @@
     thermo->correct();
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+    dimensionedScalar totalVolume = sum(mesh.V());
+

applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H

@@ -3,6 +3,10 @@
 
     rho = thermo->rho();
 
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo->rho() -= psi*p;
+
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
 
@@ -23,7 +27,7 @@
     {
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phi)
           - fvm::laplacian(rhorUAf, p)
         );
@@ -43,6 +47,9 @@
         }
     }
 
+    // Second part of thermodynamic density update
+    thermo->rho() += psi*p;
+
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
@@ -55,8 +62,10 @@
     // to obey overall mass continuity
     if (closedVolume)
     {
-        p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-            /fvc::domainIntegrate(thermo->psi());
-        rho = thermo->rho();
+        p +=
+            (initialMass - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
+        thermo->rho() = psi*p;
+        rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume;
     }
 }