chtMultiRegionFoam: updated thermodynamics

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"
@@ -121,7 +121,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info << "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicPsiThermo.H"
+#include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

     // 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());
@@ -28,7 +28,7 @@
         thermoFluid.set
         (
             i,
-            basicPsiThermo::New(fluidRegions[i]).ptr()
+            basicRhoThermo::New(fluidRegions[i]).ptr()
         );
 
         Info<< "    Adding to rhoFluid\n" << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

 {
-    // From buoyantSimpleFoam
     rho = thermo.rho();
     rho = max(rho, rhoMin[i]);
     rho = min(rho, rhoMax[i]);
@@ -13,6 +12,8 @@
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p_rgh);
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
     phi -= buoyancyPhi;
@@ -25,7 +26,11 @@
             fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
 
         p_rghEqn.solve();
 
@@ -50,10 +55,10 @@
 
     // For closed-volume cases adjust the pressure level
     // 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;
         p_rgh = p - rho*gh;
     }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

     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];

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C

@@ -34,16 +34,13 @@ Foam::scalar Foam::compressibleCourantNo
     const surfaceScalarField& phi
 )
 {
-    scalar CoNum = 0.0;
-    scalar meanCoNum = 0.0;
-
     scalarField sumPhi =
         fvc::surfaceSum(mag(phi))().internalField()
        /rho.internalField();
 
-    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+    scalar CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 
-    meanCoNum =
+    scalar meanCoNum =
         0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
 
     Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

     // 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;

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

 {
     bool closedVolume = p_rgh.needReference();
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     rho = thermo.rho();
 
@@ -19,34 +21,48 @@
 
     phi = phiU - rhorAUf*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(rhorAUf, 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(rhorAUf, 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;
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

     fvMesh& mesh = fluidRegions[i];
 
-    basicPsiThermo& thermo = thermoFluid[i];
+    basicRhoThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H

-    const dictionary& piso = solidRegions[i].solutionDict().subDict("PISO");
-
-    const int nNonOrthCorr =
-        piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -103,7 +103,7 @@
         fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
     }
 
-    Info<< "Flow time scale min/max = "
+    Info<< "Smoothed flow time scale min/max = "
         << gMin(1/rDeltaT.internalField())
         << ", " << gMax(1/rDeltaT.internalField()) << endl;
 
@@ -116,13 +116,12 @@
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
-        Info<< "Damping rDeltaT" << endl;
         rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff);
-    }
 
-    Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.internalField())
-        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.internalField())
+            << ", " << gMax(1/rDeltaT.internalField()) << endl;
+    }
 
     label nAlphaSubCycles
     (

src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C

@@ -32,6 +32,7 @@ Description
 #include "makeBasicMixture.H"
 
 #include "perfectGas.H"
+#include "incompressible.H"
 
 #include "eConstThermo.H"
 
@@ -42,6 +43,10 @@ Description
 #include "constTransport.H"
 #include "sutherlandTransport.H"
 
+#include "icoPolynomial.H"
+#include "hPolynomialThermo.H"
+#include "polynomialTransport.H"
+
 #include "pureMixture.H"
 
 #include "addToRunTimeSelectionTable.H"
@@ -94,6 +99,20 @@ makeBasicMixture
     perfectGas
 );
 
+makeBasicMixture
+(
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    incompressible
+);
+
+makeBasicPolyMixture
+(
+    pureMixture,
+    3
+);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/basic/mixtures/basicMixture/makeBasicMixture.H

@@ -44,6 +44,30 @@ defineTemplateTypeNameAndDebugWithName                                        \
     (Mixture##Transport##Thermo##EqnOfState,                                  \
     #Mixture"<"#Transport"<specieThermo<"#Thermo"<"#EqnOfState">>>>", 0)
 
+
+#define makeBasicPolyMixture(Mixture,Order)                                   \
+                                                                              \
+typedef polynomialTransport                                                   \
+<                                                                             \
+    specieThermo                                                              \
+    <                                                                         \
+        hPolynomialThermo                                                     \
+        <                                                                     \
+            icoPolynomial<Order>,                                             \
+            Order                                                             \
+            >                                                                 \
+    >,                                                                        \
+    Order                                                                     \
+> icoPoly##Order##ThermoPhysics;                                              \
+                                                                              \
+typedef Mixture<icoPoly##Order##ThermoPhysics>                                \
+    Mixture##icoPoly##Order##ThermoPhysics;                                   \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+    (Mixture##icoPoly##Order##ThermoPhysics,                                  \
+    #Mixture"<icoPoly"#Order"ThermoPhysics>", 0)
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/thermophysicalModels/basic/rhoThermo/basicRhoThermo/makeBasicRhoThermo.H

@@ -58,6 +58,39 @@ addToRunTimeSelectionTable                                                    \
 )
 
 
+#define makeBasicRhoPolyThermo(Cthermo,Mixture,Order)                         \
+                                                                              \
+typedef polynomialTransport                                                   \
+<                                                                             \
+    specieThermo                                                              \
+    <                                                                         \
+        hPolynomialThermo                                                     \
+        <                                                                     \
+            icoPolynomial<Order>,                                             \
+            Order                                                             \
+            >                                                                 \
+    >,                                                                        \
+    Order                                                                     \
+> icoPoly##Order##ThermoPhysics;                                              \
+                                                                              \
+typedef Cthermo<Mixture<icoPoly##Order##ThermoPhysics> >                      \
+    Cthermo##Mixture##icoPoly##Order##ThermoPhysics;                          \
+                                                                              \
+defineTemplateTypeNameAndDebugWithName                                        \
+(                                                                             \
+    Cthermo##Mixture##icoPoly##Order##ThermoPhysics,                          \
+    #Cthermo"<"#Mixture"<icoPoly"#Order"ThermoPhysics>>",                     \
+    0                                                                         \
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicRhoThermo,                                                           \
+    Cthermo##Mixture##icoPoly##Order##ThermoPhysics,                          \
+    fvMesh                                                                    \
+)
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/thermophysicalModels/basic/rhoThermo/hRhoThermo/hRhoThermos.C

@@ -26,6 +26,7 @@ License
 #include "makeBasicRhoThermo.H"
 
 #include "perfectGas.H"
+#include "incompressible.H"
 
 #include "hConstThermo.H"
 #include "janafThermo.H"
@@ -34,6 +35,10 @@ License
 #include "constTransport.H"
 #include "sutherlandTransport.H"
 
+#include "icoPolynomial.H"
+#include "hPolynomialThermo.H"
+#include "polynomialTransport.H"
+
 #include "hRhoThermo.H"
 #include "pureMixture.H"
 
@@ -71,6 +76,21 @@ makeBasicRhoThermo
     perfectGas
 );
 
+makeBasicRhoThermo
+(
+    hRhoThermo,
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    incompressible
+);
+
+makeBasicRhoPolyThermo
+(
+    hRhoThermo,
+    pureMixture,
+    3
+);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/specie/Make/files

@@ -6,6 +6,7 @@ reactions = reaction/reactions
 $(atomicWeights)/atomicWeights.C
 $(specie)/specie.C
 $(equationOfState)/perfectGas/perfectGas.C
+$(equationOfState)/incompressible/incompressible.C
 $(reactions)/makeChemkinReactions.C
 $(reactions)/makeReactionThermoReactions.C
 $(reactions)/makeLangmuirHinshelwoodReactions.C

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/T

@@ -23,6 +23,7 @@ boundaryField
     ".*"
     {
         type            calculated;
+        value           uniform 300;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/U

@@ -23,6 +23,7 @@ boundaryField
     ".*"
     {
         type            calculated;
+        value           uniform (0.01 0 0);
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/epsilon

@@ -23,6 +23,7 @@ boundaryField
     ".*"
     {
         type            calculated;
+        value           uniform 0.01;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/k

@@ -23,6 +23,7 @@ boundaryField
     ".*"
     {
         type            calculated;
+        value           uniform 0.1;
     }
 }