chtMultiRegionFoam: Changed thermodynamics from psi-based to rho-base…

…d to support liquids
OpenCFD · Oct 14, 2010 · 16a7cff · 16a7cff
1 parent b32f406
commit 16a7cff
Show file tree

Hide file tree

Showing 4 changed files with 40 additions and 24 deletions.
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicPsiThermo.H"
+#include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
@@ -118,7 +118,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info << "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -1,5 +1,5 @@
     // Initialise fluid field pointer lists
-    PtrList<basicPsiThermo> thermoFluid(fluidRegions.size());
+    PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
     PtrList<volScalarField> rhoFluid(fluidRegions.size());
     PtrList<volScalarField> KFluid(fluidRegions.size());
     PtrList<volVectorField> UFluid(fluidRegions.size());
@@ -23,7 +23,7 @@
         thermoFluid.set
         (
             i,
-            basicPsiThermo::New(fluidRegions[i]).ptr()
+            basicRhoThermo::New(fluidRegions[i]).ptr()
         );
 
         Info<< "    Adding to rhoFluid\n" << endl;

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -1,5 +1,7 @@
 {
     bool closedVolume = p_rgh.needReference();
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     rho = thermo.rho();
 
@@ -19,34 +21,48 @@
 
     phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix p_rghEqn
+        fvScalarMatrix p_rghDDtEqn
         (
-            fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p_rgh)
         );
 
-        p_rghEqn.solve
-        (
-            mesh.solver
+        // Thermodynamic density needs to be updated by psi*d(p) after the
+        // pressure solution - done in 2 parts. Part 1:
+        thermo.rho() -= psi*p_rgh;
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix p_rghEqn
+            (
+                p_rghDDtEqn
+              - fvm::laplacian(rhorUAf, p_rgh)
+            );
+
+            p_rghEqn.solve
             (
-                p_rgh.select
+                mesh.solver
                 (
+                    p_rgh.select
                     (
-                        oCorr == nOuterCorr-1
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
+                        (
+                            oCorr == nOuterCorr-1
+                         && corr == nCorr-1
+                         && nonOrth == nNonOrthCorr
+                        )
                     )
                 )
-            )
-        );
+            );
 
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi += p_rghEqn.flux();
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += p_rghEqn.flux();
+            }
         }
+
+        // Second part of thermodynamic density update
+        thermo.rho() += psi*p_rgh;
     }
 
     // Correct velocity field
@@ -66,10 +82,10 @@
 
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
-    if (closedVolume)
+    if (closedVolume && compressible)
     {
-        p += (initialMass - fvc::domainIntegrate(psi*p))
-            /fvc::domainIntegrate(psi);
+        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+            /compressibility;
         rho = thermo.rho();
         p_rgh = p - rho*gh;
     }

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -1,6 +1,6 @@
     const fvMesh& mesh = fluidRegions[i];
 
-    basicPsiThermo& thermo = thermoFluid[i];
+    basicRhoThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];