incompressible/simpleSRFFoam, solver + tutorials

applications/solvers/incompressible/simpleSRFFoam/UEqn.H

@@ -0,0 +1,14 @@
+    // Solve the momentum equation
+
+    tmp<fvVectorMatrix> HUrelEqn
+    (
+        fvm::div(phi, Urel)
+      + turbulence->divDevReff(Urel)
+      + SRF->Su()
+    );
+
+    // Get under-relaxation factor
+    const scalar UUrf = mesh.solutionDict().relaxationFactor(Urel.name());
+
+    // Momentum solution
+    solve(relax(HUrelEqn(), UUrf) == -fvc::grad(p));

applications/solvers/incompressible/simpleSRFFoam/createFields.H

@@ -71,4 +71,3 @@
         ),
         Urel + SRF->U()
     );
-

applications/solvers/incompressible/simpleSRFFoam/pEqn.H

@@ -0,0 +1,48 @@
+    p.boundaryField().updateCoeffs();
+
+    // Prepare clean 1/Ap without contribution from under-relaxation
+    // HJ, 26/Oct/2015
+    volScalarField rUA
+    (
+        "(1|A(Urel))",
+        1/HUrelEqn().A()
+    );
+
+    // Store velocity under-relaxation point before using U for the flux
+    // precursor
+    Urel.storePrevIter();
+
+    Urel = rUA*HUrelEqn().H();
+    HUrelEqn.clear();
+    phi = fvc::interpolate(Urel) & mesh.Sf();
+
+    // Global flux balance
+    adjustPhi(phi, Urel, p);
+
+    // Non-orthogonal pressure corrector loop
+    for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rUA, p) == fvc::div(phi)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= pEqn.flux();
+        }
+    }
+
+#   include "continuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    // Momentum corrector
+    // Note: since under-relaxation does not change aU, H/a in U can be
+    // re-used.  HJ, 22/Jan/2016
+    Urel = UUrf*(Urel - rUA*fvc::grad(p)) + (1 - UUrf)*Urel.prevIter();
+    Urel.correctBoundaryConditions();

applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C

@@ -25,8 +25,12 @@ Application
     simpleSRFFoam
 
 Description
-    Steady-state solver for incompressible, turbulent flow of non-Newtonian
+    Steady-state solver for incompressible, turbulent flow of Newtonian
     fluids with single rotating frame.
+    Consistent formulation without time-step and relaxation dependence by Jasak
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
 
 \*---------------------------------------------------------------------------*/
 
@@ -47,8 +51,6 @@ int main(int argc, char *argv[])
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 
-    //mesh.clearPrimitives();
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -63,50 +65,8 @@ int main(int argc, char *argv[])
 
         // Pressure-velocity SIMPLE corrector
         {
-            // Momentum predictor
-            tmp<fvVectorMatrix> UrelEqn
-            (
-                fvm::div(phi, Urel)
-              + turbulence->divDevReff(Urel)
-              + SRF->Su()
-            );
-
-            UrelEqn().relax();
-
-            solve(UrelEqn() == -fvc::grad(p));
-
-            p.boundaryField().updateCoeffs();
-            volScalarField AUrel = UrelEqn().A();
-            Urel = UrelEqn().H()/AUrel;
-            UrelEqn.clear();
-            phi = fvc::interpolate(Urel) & mesh.Sf();
-            adjustPhi(phi, Urel, p);
-
-            // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::laplacian(1.0/AUrel, p) == fvc::div(phi)
-                );
-
-                pEqn.setReference(pRefCell, pRefValue);
-                pEqn.solve();
-
-                if (nonOrth == nNonOrthCorr)
-                {
-                    phi -= pEqn.flux();
-                }
-            }
-
-#           include "continuityErrs.H"
-
-            // Explicitly relax pressure for momentum corrector
-            p.relax();
-
-            // Momentum corrector
-            Urel -= fvc::grad(p)/AUrel;
-            Urel.correctBoundaryConditions();
+#           include "UEqn.H"
+#           include "pEqn.H"
         }
 
         turbulence->correct();

tutorials/incompressible/simpleSRFFoam/mixer/system/fvSolution

@@ -25,8 +25,8 @@ solvers
     };
     Urel
     {
-        solver           PBiCG;
-        preconditioner   DILU;
+        solver           smoothSolver;
+        smoother         GaussSeidel;
         tolerance        1e-05;
         relTol           0.1;
     };
@@ -75,7 +75,7 @@ SIMPLE
 relaxationFactors
 {
     p               0.3;
-    Urel            0.7;
+    Urel            0.5;
     k               0.7;
     epsilon         0.7;
     omega           0.7;