improved flux calcs + more clean-up

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

@@ -10,5 +10,16 @@
 
     if (momentumPredictor)
     {
-        solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh);
+         solve
+         (
+            UEqn()
+         ==
+           -fvc::reconstruct
+            (
+                (
+                    fvc::snGrad(pd)
+                  + ghf*fvc::snGrad(rho)
+                ) * mesh.magSf()
+            )
+        );
     }

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

@@ -7,8 +7,8 @@
             (
                 fluidRegions[regionI],
                 runTime,
-                rhof[regionI],
-                phif[regionI]
+                rhoFluid[regionI],
+                phiFluid[regionI]
             ),
             CoNum
         );

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

@@ -1,15 +1,16 @@
     // Initialise fluid field pointer lists
-    PtrList<basicThermo> thermof(fluidRegions.size());
-    PtrList<volScalarField> rhof(fluidRegions.size());
-    PtrList<volScalarField> Kf(fluidRegions.size());
-    PtrList<volVectorField> Uf(fluidRegions.size());
-    PtrList<surfaceScalarField> phif(fluidRegions.size());
+    PtrList<basicThermo> thermoFluid(fluidRegions.size());
+    PtrList<volScalarField> rhoFluid(fluidRegions.size());
+    PtrList<volScalarField> KFluid(fluidRegions.size());
+    PtrList<volVectorField> UFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-    PtrList<volScalarField> DpDtf(fluidRegions.size());
-    PtrList<volScalarField> ghf(fluidRegions.size());
-    PtrList<volScalarField> pdf(fluidRegions.size());
+    PtrList<volScalarField> DpDtFluid(fluidRegions.size());
+    PtrList<volScalarField> ghFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
+    PtrList<volScalarField> pdFluid(fluidRegions.size());
 
-    List<scalar> initialMassf(fluidRegions.size());
+    List<scalar> initialMassFluid(fluidRegions.size());
 
     dimensionedScalar pRef
     (
@@ -24,8 +25,8 @@
         Info<< "*** Reading fluid mesh thermophysical properties for region "
             << fluidRegions[i].name() << nl << endl;
 
-        Info<< "    Adding to pdf\n" << endl;
-        pdf.set
+        Info<< "    Adding to pdFluid\n" << endl;
+        pdFluid.set
         (
             i,
             new volScalarField
@@ -42,16 +43,15 @@
             )
         );
 
-        Info<< "    Adding to thermof\n" << endl;
-
-        thermof.set
+        Info<< "    Adding to thermoFluid\n" << endl;
+        thermoFluid.set
         (
             i,
             basicThermo::New(fluidRegions[i]).ptr()
         );
 
-        Info<< "    Adding to rhof\n" << endl;
-        rhof.set
+        Info<< "    Adding to rhoFluid\n" << endl;
+        rhoFluid.set
         (
             i,
             new volScalarField
@@ -64,12 +64,12 @@
                     IOobject::NO_READ,
                     IOobject::AUTO_WRITE
                 ),
-                thermof[i].rho()
+                thermoFluid[i].rho()
             )
         );
 
-        Info<< "    Adding to Kf\n" << endl;
-        Kf.set
+        Info<< "    Adding to KFluid\n" << endl;
+        KFluid.set
         (
             i,
             new volScalarField
@@ -82,12 +82,12 @@
                     IOobject::NO_READ,
                     IOobject::NO_WRITE
                 ),
-                thermof[i].Cp()*thermof[i].alpha()
+                thermoFluid[i].Cp()*thermoFluid[i].alpha()
             )
         );
 
-        Info<< "    Adding to Uf\n" << endl;
-        Uf.set
+        Info<< "    Adding to UFluid\n" << endl;
+        UFluid.set
         (
             i,
             new volVectorField
@@ -104,8 +104,8 @@
             )
         );
 
-        Info<< "    Adding to phif\n" << endl;
-        phif.set
+        Info<< "    Adding to phiFluid\n" << endl;
+        phiFluid.set
         (
             i,
             new surfaceScalarField
@@ -118,7 +118,7 @@
                     IOobject::READ_IF_PRESENT,
                     IOobject::AUTO_WRITE
                 ),
-                linearInterpolate(rhof[i]*Uf[i])
+                linearInterpolate(rhoFluid[i]*UFluid[i])
                     & fluidRegions[i].Sf()
             )
         );
@@ -131,16 +131,16 @@
             (
                 compressible::turbulenceModel::New
                 (
-                    rhof[i],
-                    Uf[i],
-                    phif[i],
-                    thermof[i]
+                    rhoFluid[i],
+                    UFluid[i],
+                    phiFluid[i],
+                    thermoFluid[i]
                 )
             ).ptr()
         );
 
-        Info<< "    Adding to DpDtf\n" << endl;
-        DpDtf.set
+        Info<< "    Adding to DpDtFluid\n" << endl;
+        DpDtFluid.set
         (
             i,
             new volScalarField
@@ -150,9 +150,9 @@
                     surfaceScalarField
                     (
                         "phiU",
-                        phif[i]/fvc::interpolate(rhof[i])
+                        phiFluid[i]/fvc::interpolate(rhoFluid[i])
                     ),
-                    thermof[i].p()
+                    thermoFluid[i].p()
                 )
             )
         );
@@ -162,8 +162,8 @@
                 ("environmentalProperties");
         dimensionedVector g(environmentalProperties.lookup("g"));
 
-        Info<< "    Adding to ghf\n" << endl;
-        ghf.set
+        Info<< "    Adding to ghFluid\n" << endl;
+        ghFluid.set
         (
             i,
             new volScalarField
@@ -172,12 +172,21 @@
                  g & fluidRegions[i].C()
             )
         );
+        ghfFluid.set
+        (
+            i,
+            new surfaceScalarField
+            (
+                "ghf",
+                 g & fluidRegions[i].Cf()
+            )
+        );
 
         Info<< "    Updating p from pd\n" << endl;
-        thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef;
-        thermof[i].correct();
+        thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef;
+        thermoFluid[i].correct();
 
-        initialMassf[i] = fvc::domainIntegrate(rhof[i]).value();
+        initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
     }
 
 

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

@@ -1,24 +1,31 @@
 {
-    bool closedVolume = false;
+    bool closedVolume = pd.needReference();
 
     rho = thermo.rho();
 
     volScalarField rUA = 1.0/UEqn().A();
+    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+
     U = rUA*UEqn().H();
 
-    phi =
+    surfaceScalarField phiU
+    (
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         )
-      - fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf();
+    );
+
+    phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pdEqn
         (
-            fvc::ddt(rho)
+            fvm::ddt(psi, pd)
+          + fvc::ddt(psi)*pRef
+          + fvc::ddt(psi, rho)*gh
           + fvc::div(phi)
           - fvm::laplacian(rho*rUA, pd)
         );
@@ -38,8 +45,12 @@
         }
     }
 
+    // Correct velocity field
+    U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
+    U.correctBoundaryConditions();
+
     // Update pressure field (including bc)
-    thermo.p() == pd + rho*gh + pRef;
+    p == pd + rho*gh + pRef;
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     // Solve continuity
@@ -48,23 +59,17 @@
     // Update continuity errors
 #   include "compressibleContinuityErrors.H"
 
-    // Correct velocity field
-    U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
-    U.correctBoundaryConditions();
-
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
     {
-        thermo.p() +=
-            (
-                dimensionedScalar("massIni", dimMass, initialMassf[i])
-              - fvc::domainIntegrate(thermo.psi()*thermo.p())
-            )/fvc::domainIntegrate(thermo.psi());
-        pd == thermo.p() - (rho*gh + pRef);
+        p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi);
         rho = thermo.rho();
     }
 
     // Update thermal conductivity
-    K = thermo.Cp()*turb.alphaEff();
+    K = thermoFluid[i].Cp()*turb.alphaEff();
+
+    // Update pd (including bc)
+    pd == p - (rho*gh + pRef);
 }