From 56a81f1825384f96737141118e65ca906e07e9fd Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 26 Sep 2019 11:21:23 +0100
Subject: [PATCH] Updated version of buoyantBoussinesqSimpleFoam

---
 .../buoyantBoussinesqSimpleFoam/TEqn.H        |  3 +-
 .../buoyantBoussinesqSimpleFoam/UEqn.H        | 13 ++----
 .../buoyantBoussinesqSimpleFoam.C             | 18 ++++----
 .../buoyantBoussinesqSimpleFoam/pEqn.H        | 42 ++++++++++++-------
 4 files changed, 40 insertions(+), 36 deletions(-)

diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
index 3bdeea045..73887cb96 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
@@ -8,12 +8,11 @@
     fvScalarMatrix TEqn
     (
         fvm::div(phi, T)
-      - fvm::Sp(fvc::div(phi), T)
+      + fvm::SuSp(-fvc::div(phi), T)
       - fvm::laplacian(kappaEff, T)
     );
 
     TEqn.relax();
-
     TEqn.solve();
 
     rhok = 1.0 - beta*(T - TRef);
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index 4fcf603b1..3155621db 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -6,17 +6,10 @@
       + turbulence->divDevReff()
     );
 
-    UEqn().relax();
-
     solve
     (
-        UEqn()
+        relax(UEqn())
       ==
-        fvc::reconstruct
-        (
-            (
-                fvc::interpolate(rhok)*(g & mesh.Sf())
-              - fvc::snGrad(p)*mesh.magSf()
-            )
-        )
+       -fvc::grad(p)
+      + fvc::reconstruct(fvc::interpolate(rhok)*(g & mesh.Sf()))
     );
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index e49c7a40b..560f9fcbb 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -54,15 +54,15 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMesh.H"
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
 
     simpleControl simple(mesh);
 
-    #include "readGravitationalAcceleration.H"
-    #include "createFields.H"
-    #include "initContinuityErrs.H"
+#   include "readGravitationalAcceleration.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,9 +74,9 @@ int main(int argc, char *argv[])
 
         // Pressure-velocity SIMPLE corrector
         {
-            #include "UEqn.H"
-            #include "TEqn.H"
-            #include "pEqn.H"
+#           include "TEqn.H"
+#           include "UEqn.H"
+#           include "pEqn.H"
         }
 
         turbulence->correct();
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 256a336f3..c3d59b0d6 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -1,22 +1,31 @@
 {
-    volScalarField rUA("rUA", 1.0/UEqn().A());
-    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+    p.boundaryField().updateCoeffs();
 
-    U = rUA*UEqn().H();
+    volScalarField rAU("rAU", 1.0/UEqn().A());
+    surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
+
+    U = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(U) & mesh.Sf();
+    simple.calcSteadyConsistentFlux(phi, U);
     adjustPhi(phi, U, p);
 
     surfaceScalarField buoyancyPhi =
-        rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
+        rAUf*fvc::interpolate(rhok)*(g & mesh.Sf());
     phi += buoyancyPhi;
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUAf, p) == fvc::div(phi)
+            fvm::laplacian
+            (
+                rAUf/simple.aCoeff(U.name()),
+                p,
+                "laplacian(rAU," + p.name() + ')'
+            )
+         ==
+            fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -27,16 +36,19 @@
         {
             // Calculate the conservative fluxes
             phi -= pEqn.flux();
-
-            // Explicitly relax pressure for momentum corrector
-            p.relax();
-
-            // Correct the momentum source with the pressure gradient flux
-            // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf);
-            U.correctBoundaryConditions();
         }
+
+        // U += rAU*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rAUf);
+        // U.correctBoundaryConditions();
     }
 
-    #include "continuityErrs.H"
+#   include "continuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    // Correct the momentum source with the pressure gradient flux
+    // calculated from the relaxed pressure
+    U += rAU*(fvc::reconstruct(buoyancyPhi/rAUf) - fvc::grad(p));
+    U.correctBoundaryConditions();
 }