From 1b6f46fae03e985fd983b7dd747287ee1d6db37e Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Fri, 29 May 2015 11:30:40 +0100
Subject: [PATCH] simpleFoam/SRFSimpleFoam: Added support for SIMPLEC SIMPLEC
 (SIMPLE-consistent) is selected by setting "consistent" option true/yes:

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent yes;
}

which relaxes the pressure in a "consistent" manner and additional
relaxation of the pressure is not generally necessary.  In addition
convergence of the p-U system is better and reliable with less
aggressive relaxation of the momentum equation, e.g. for the motorbike
tutorial:

relaxationFactors
{
    equations
    {
        U               0.9;
        k               0.7;
        omega           0.7;
    }
}

The cost per iteration is marginally higher but the convergence rate is
better so the number of iterations can be reduced.

The SIMPLEC algorithm also provides benefit for cases with large
body-forces, e.g. SRF, see tutorials/incompressible/SRFSimpleFoam/mixer
and feature request http://www.openfoam.org/mantisbt/view.php?id=1714
---
 .../simpleFoam/SRFSimpleFoam/pEqn.H           | 17 +++++++++--
 .../solvers/incompressible/simpleFoam/pEqn.H  | 19 ++++++++----
 .../solutionControl/solutionControl.C         | 18 +++++++----
 .../solutionControl/solutionControl.H         | 10 ++++++-
 .../solutionControl/solutionControlI.H        |  8 ++++-
 .../SRFSimpleFoam/mixer/system/fvSolution     | 12 ++------
 .../SRFSimpleFoam/mixer/turbulenceProperties  | 30 -------------------
 .../simpleFoam/motorBike/system/fvSolution    |  7 ++---
 8 files changed, 60 insertions(+), 61 deletions(-)
 delete mode 100644 tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties

diff --git a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
index b9fa52c0f..50691c015 100644
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@@ -2,17 +2,28 @@
     volScalarField rAUrel(1.0/UrelEqn().A());
     volVectorField HbyA("HbyA", Urel);
     HbyA = rAUrel*UrelEqn().H();
-    UrelEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
     adjustPhi(phiHbyA, Urel, p);
 
+    tmp<volScalarField> rAtUrel(rAUrel);
+
+    if (simple.consistent())
+    {
+        rAtUrel = 1.0/(1.0/rAUrel - UrelEqn().H1());
+        phiHbyA +=
+            fvc::interpolate(rAtUrel() - rAUrel)*fvc::snGrad(p)*mesh.magSf();
+        HbyA -= (rAUrel - rAtUrel())*fvc::grad(p);
+    }
+
+    UrelEqn.clear();
+
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAUrel, p) == fvc::div(phiHbyA)
+            fvm::laplacian(rAtUrel(), p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -31,7 +42,7 @@
     p.relax();
 
     // Momentum corrector
-    Urel = HbyA - rAUrel*fvc::grad(p);
+    Urel = HbyA - rAtUrel()*fvc::grad(p);
     Urel.correctBoundaryConditions();
     fvOptions.correct(Urel);
 }
diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H
index 4d669c8db..fddb3b05d 100644
--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@@ -2,20 +2,29 @@
     volScalarField rAU(1.0/UEqn().A());
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
-    UEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
-
     fvOptions.makeRelative(phiHbyA);
-
     adjustPhi(phiHbyA, U, p);
 
+    tmp<volScalarField> rAtU(rAU);
+
+    if (simple.consistent())
+    {
+        rAtU = 1.0/(1.0/rAU - UEqn().H1());
+        phiHbyA +=
+            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
+        HbyA -= (rAU - rAtU())*fvc::grad(p);
+    }
+
+    UEqn.clear();
+
     // Non-orthogonal pressure corrector loop
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -34,7 +43,7 @@
     p.relax();
 
     // Momentum corrector
-    U = HbyA - rAU*fvc::grad(p);
+    U = HbyA - rAtU()*fvc::grad(p);
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 }
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
index 8448153a6..8e3dd8907 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,16 +37,21 @@ namespace Foam
 
 void Foam::solutionControl::read(const bool absTolOnly)
 {
-    const dictionary& solnDict = this->dict();
+    const dictionary& solutionDict = this->dict();
 
     // Read solution controls
     nNonOrthCorr_ =
-        solnDict.lookupOrDefault<label>("nNonOrthogonalCorrectors", 0);
-    momentumPredictor_ = solnDict.lookupOrDefault("momentumPredictor", true);
-    transonic_ = solnDict.lookupOrDefault("transonic", false);
+        solutionDict.lookupOrDefault<label>("nNonOrthogonalCorrectors", 0);
+    momentumPredictor_ =
+        solutionDict.lookupOrDefault("momentumPredictor", true);
+    transonic_ = solutionDict.lookupOrDefault("transonic", false);
+    consistent_ = solutionDict.lookupOrDefault("consistent", false);
 
     // Read residual information
-    const dictionary residualDict(solnDict.subOrEmptyDict("residualControl"));
+    const dictionary residualDict
+    (
+        solutionDict.subOrEmptyDict("residualControl")
+    );
 
     DynamicList<fieldData> data(residualControl_);
 
@@ -177,6 +182,7 @@ Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
     nNonOrthCorr_(0),
     momentumPredictor_(true),
     transonic_(false),
+    consistent_(false),
     corr_(0),
     corrNonOrtho_(0)
 {}
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
index e701813d4..9cd6cbbf1 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,6 +83,10 @@ protected:
             //- Flag to indicate to solve using transonic algorithm
             bool transonic_;
 
+            //- Flag to indicate to relax pressure using the
+            //  "consistent" approach of SIMPLEC
+            bool consistent_;
+
 
         // Evolution
 
@@ -169,6 +173,10 @@ public:
             //- Flag to indicate to solve using transonic algorithm
             inline bool transonic() const;
 
+            //- Flag to indicate to relax pressure using the
+            //  "consistent" approach of SIMPLEC
+            inline bool consistent() const;
+
 
         // Evolution
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
index b9d975cfb..61a647e40 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,6 +67,12 @@ inline bool Foam::solutionControl::transonic() const
 }
 
 
+inline bool Foam::solutionControl::consistent() const
+{
+    return consistent_;
+}
+
+
 inline bool Foam::solutionControl::correctNonOrthogonal()
 {
     corrNonOrtho_++;
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution b/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution
index 99306fab0..e092956f3 100644
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution
+++ b/tutorials/incompressible/SRFSimpleFoam/mixer/system/fvSolution
@@ -37,22 +37,14 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
+    consistent          yes;
 }
 
 relaxationFactors
 {
-    fields
-    {
-        p               0.3;
-    }
     equations
     {
-        Urel            0.7;
-        k               0.7;
-        epsilon         0.7;
-        omega           0.7;
-        R               0.7;
-        nuTilda         0.7;
+        ".*"            0.9;
     }
 }
 
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties b/tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties
deleted file mode 100644
index 4406361c6..000000000
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/turbulenceProperties
+++ /dev/null
@@ -1,30 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      turbulenceProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-simulationType RAS;
-
-RAS
-{
-    RASModel        kOmegaSST;
-
-    turbulence      on;
-
-    printCoeffs     on;
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
index 4f0f515ac..6f9e2f67b 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/fvSolution
@@ -66,6 +66,7 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
+    consistent yes;
 }
 
 potentialFlow
@@ -75,13 +76,9 @@ potentialFlow
 
 relaxationFactors
 {
-    fields
-    {
-        p               0.3;
-    }
     equations
     {
-        U               0.7;
+        U               0.9;
         k               0.7;
         omega           0.7;
     }