applying improved flux-velocity correspondence

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -1,4 +1,4 @@
-    // Solve the Momentum equation
+    // Solve the momentum equation
 
     tmp<fvVectorMatrix> UEqn
     (
@@ -13,8 +13,13 @@
     (
         UEqn()
       ==
-       -fvc::grad(pd)
-      + beta*gh*fvc::grad(T)
+       -fvc::reconstruct
+        (
+            (
+                fvc::snGrad(pd)
+              - betaghf*fvc::snGrad(T)
+            ) * mesh.magSf()
+        )
     ).initialResidual();
 
     maxResidual = max(eqnResidual, maxResidual);

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H

@@ -53,8 +53,8 @@
         incompressible::RASModel::New(U, phi, laminarTransport)
     );
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
+    Info<< "Calculating field beta*(g.h)\n" << endl;
+    surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
 
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H

@@ -1,41 +1,49 @@
-volScalarField rUA = 1.0/UEqn().A();
-U = rUA*UEqn().H();
-UEqn.clear();
-
-phi = fvc::interpolate(U) & mesh.Sf();
-adjustPhi(phi, U, pd);
-phi += fvc::interpolate(beta*gh*rUA)*fvc::snGrad(T)*mesh.magSf();
-
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-    fvScalarMatrix pdEqn
-    (
-        fvm::laplacian(rUA, pd) == fvc::div(phi)
-    );
+    volScalarField rUA("rUA", 1.0/UEqn().A());
+    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
 
-    pdEqn.setReference(pdRefCell, pdRefValue);
+    U = rUA*UEqn().H();
+    UEqn.clear();
 
-    // retain the residual from the first iteration
-    if (nonOrth == 0)
+    phi = fvc::interpolate(U) & mesh.Sf();
+    adjustPhi(phi, U, pd);
+    surfaceScalarField buoyancyPhi = -betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
+    phi -= buoyancyPhi;
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        eqnResidual = pdEqn.solve().initialResidual();
-        maxResidual = max(eqnResidual, maxResidual);
-    }
-    else
-    {
-        pdEqn.solve();
+        fvScalarMatrix pdEqn
+        (
+            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+        );
+
+        pdEqn.setReference(pdRefCell, pdRefValue);
+
+        // retain the residual from the first iteration
+        if (nonOrth == 0)
+        {
+            eqnResidual = pdEqn.solve().initialResidual();
+            maxResidual = max(eqnResidual, maxResidual);
+        }
+        else
+        {
+            pdEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            // Calculate the conservative fluxes
+            phi -= pdEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            pd.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rUAf);
+            U.correctBoundaryConditions();
+        }
     }
 
-    if (nonOrth == nNonOrthCorr)
-    {
-        phi -= pdEqn.flux();
-    }
+    #include "continuityErrs.H"
 }
-
-#include "continuityErrs.H"
-
-// Explicitly relax pressure for momentum corrector
-pd.relax();
-
-U -= rUA*(fvc::grad(pd) - beta*gh*fvc::grad(T));
-U.correctBoundaryConditions();

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H

@@ -23,7 +23,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        pd + rhoEff*gh + pRef
+        pd + rhoEff*(g & mesh.C()) + pRef
     );
     p.write();
 }