ENH: Added relaxed nonOrthogonal laplacian corrector scheme

src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.C

@@ -190,6 +190,15 @@ realizableKE<BasicTurbulenceModel>::realizableKE
             1.2
         )
     ),
+    implicitKEpsilonSource_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "implicitKEpsilonSource",
+            this->coeffDict_,
+            false
+        )
+    ),
 
     k_
     (
@@ -237,6 +246,11 @@ bool realizableKE<BasicTurbulenceModel>::read()
         C2_.readIfPresent(this->coeffDict());
         sigmak_.readIfPresent(this->coeffDict());
         sigmaEps_.readIfPresent(this->coeffDict());
+        implicitKEpsilonSource_.readIfPresent
+        (
+            "implicitKEpsilonSource",
+            this->coeffDict()
+        );
 
         return true;
     }
@@ -288,6 +302,41 @@ void realizableKE<BasicTurbulenceModel>::correct()
         )
     );
 
+    tmp<fvScalarMatrix> tEpsilonSourceEq
+    (
+        new fvScalarMatrix
+        (
+            epsilon_,
+            dimVolume*rho.dimensions()*epsilon_.dimensions()/dimTime
+        )
+    );
+
+
+    tmp<volScalarField> Cmu(rCmu(tgradU(), S2, magS));
+
+    if (implicitKEpsilonSource_)
+    {
+        tmp<volScalarField> tEpsSp
+        (
+           - C1*alpha*rho*eta/sqrt(this->nut_/Cmu())
+           + alpha*rho*C2_/(sqrt(this->nut_/Cmu()) + sqrt(this->nu()))
+        );
+
+        tEpsilonSourceEq.ref() +=
+            0.5*tEpsSp()*pow(epsilon_, 1.5)
+          - fvm::Sp(1.5*tEpsSp()*sqrt(epsilon_), epsilon_);
+    }
+    else
+    {
+        tEpsilonSourceEq.ref() +=
+            C1*alpha*rho*magS*epsilon_
+          - fvm::Sp
+            (
+                C2_*alpha*rho*epsilon_/(k_ + sqrt(nuLimited*epsilon_)),
+                epsilon_
+            );
+    }
+
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
     (
@@ -295,12 +344,7 @@ void realizableKE<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, epsilon_)
       - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
      ==
-        C1*alpha*rho*magS*epsilon_
+        tEpsilonSourceEq()
-      - fvm::Sp
-        (
-            C2_*alpha*rho*epsilon_/(k_ + sqrt(nuLimited*epsilon_)),
-            epsilon_
-        )
       + epsilonSource()
       + fvOptions(alpha, rho, epsilon_)
     );
@@ -314,6 +358,26 @@ void realizableKE<BasicTurbulenceModel>::correct()
 
 
     // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> tkSourceEq
+    (
+        new fvScalarMatrix
+        (
+            k_,
+            dimVolume*rho.dimensions()*k_.dimensions()/dimTime
+        )
+    );
+
+    if (implicitKEpsilonSource_)
+    {
+        tkSourceEq.ref() +=
+            alpha*rho*G + alpha*alpha*rho*rho*Cmu()/this->nut_*sqr(k_)
+          - fvm::Sp(2.0*alpha*alpha*rho*rho*Cmu()*k_/this->nut_, k_);
+    }
+    else
+    {
+        tkSourceEq.ref() += alpha*rho*G - fvm::Sp(alpha*rho*epsilon_/k_, k_);
+    }
+
 
     tmp<fvScalarMatrix> kEqn
     (
@@ -321,9 +385,8 @@ void realizableKE<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha*rho*G
       - fvm::SuSp(2.0/3.0*alpha*rho*divU, k_)
-      - fvm::Sp(alpha*rho*epsilon_/k_, k_)
+      + tkSourceEq()
       + kSource()
       + fvOptions(alpha, rho, k_)
     );

src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.H

@@ -97,6 +97,9 @@ protected:
             dimensionedScalar sigmak_;
             dimensionedScalar sigmaEps_;
 
+            //- Use implicit k-Epilon source
+            Switch implicitKEpsilonSource_;
+
 
         // Fields
 

src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C

@@ -0,0 +1,276 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "relaxedNonOrthoGaussLaplacianScheme.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
+(
+    const surfaceScalarField& gammaMagSf,
+    const surfaceScalarField& deltaCoeffs,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type>> tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
+        )
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
+    fvm.negSumDiag();
+
+    forAll(vf.boundaryField(), patchi)
+    {
+        const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
+        const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
+        const fvsPatchScalarField& pDeltaCoeffs =
+            deltaCoeffs.boundaryField()[patchi];
+
+        if (pvf.coupled())
+        {
+            fvm.internalCoeffs()[patchi] =
+                pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
+            fvm.boundaryCoeffs()[patchi] =
+               -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
+        }
+        else
+        {
+            fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
+            fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
+        }
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::gammaSnGradCorr
+(
+    const surfaceVectorField& SfGammaCorr,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tgammaSnGradCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "gammaSnGradCorr("+vf.name()+')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            SfGammaCorr.dimensions()
+           *vf.dimensions()*mesh.deltaCoeffs().dimensions()
+        )
+    );
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        tgammaSnGradCorr.ref().replace
+        (
+            cmpt,
+            fvc::dotInterpolate(SfGammaCorr, fvc::grad(vf.component(cmpt)))
+        );
+    }
+
+    return tgammaSnGradCorr;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
+    (
+        fvc::div(this->tsnGradScheme_().snGrad(vf)*mesh.magSf())
+    );
+
+    tLaplacian.ref().rename("laplacian(" + vf.name() + ')');
+
+    return tLaplacian;
+}
+
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType;
+
+    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
+
+    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
+    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
+    (
+        SfGamma & Sn
+    );
+    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
+
+    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected
+    (
+        SfGammaSn,
+        this->tsnGradScheme_().deltaCoeffs(vf),
+        vf
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    tmp<SType> tfaceFluxCorrection = gammaSnGradCorr(SfGammaCorr, vf);
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        tfaceFluxCorrection.ref() +=
+            SfGammaSn*this->tsnGradScheme_().correction(vf);
+    }
+
+    const word corrName(tfaceFluxCorrection().name());
+
+    tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));
+
+    const word oldName(corrName + "_0");
+    const scalar relax(vf.mesh().equationRelaxationFactor(oldName));
+
+    const objectRegistry& obr = vf.db();
+    if (obr.foundObject<SType>(oldName))
+    {
+        SType& oldCorrection = obr.lookupObjectRef<SType>(oldName);
+
+        trelaxedCorrection.ref() *= relax;
+        trelaxedCorrection.ref() += (1.0-relax)*oldCorrection;
+
+        oldCorrection = tfaceFluxCorrection;
+    }
+    else
+    {
+        SType* s = new SType(oldName, tfaceFluxCorrection);
+        s->store();
+    }
+
+    fvm.source() -=
+        mesh.V()
+       *fvc::div
+        (
+            trelaxedCorrection()
+        )().primitiveField();
+
+    if (mesh.fluxRequired(vf.name()))
+    {
+        fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
+    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
+    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
+    (
+        SfGamma & Sn
+    );
+    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
+    (
+        fvc::div
+        (
+            SfGammaSn*this->tsnGradScheme_().snGrad(vf)
+          + gammaSnGradCorr(SfGammaCorr, vf)
+        )
+    );
+
+    tLaplacian.ref().rename
+    (
+        "laplacian(" + gamma.name() + ',' + vf.name() + ')'
+    );
+
+    return tLaplacian;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.H

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::relaxedNonOrthoGaussLaplacianScheme
+
+Description
+    Basic second-order laplacian using face-gradients and Gauss' theorem.
+
+SourceFiles
+    relaxedNonOrthoGaussLaplacianScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef relaxedNonOrthoGaussLaplacianScheme_H
+#define relaxedNonOrthoGaussLaplacianScheme_H
+
+#include "laplacianScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+             Class relaxedNonOrthoGaussLaplacianScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class GType>
+class relaxedNonOrthoGaussLaplacianScheme
+:
+    public fv::laplacianScheme<Type, GType>
+{
+    // Private Member Functions
+
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
+        (
+            const surfaceVectorField& SfGammaCorr,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- No copy construct
+        relaxedNonOrthoGaussLaplacianScheme
+        (
+            const relaxedNonOrthoGaussLaplacianScheme&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const relaxedNonOrthoGaussLaplacianScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("relaxedNonOrthoGauss");
+
+
+    // Constructors
+
+        //- Construct null
+        relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh)
+        :
+            laplacianScheme<Type, GType>(mesh)
+        {}
+
+        //- Construct from Istream
+        relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh, Istream& is)
+        :
+            laplacianScheme<Type, GType>(mesh, is)
+        {}
+
+        //- Construct from mesh, interpolation and snGradScheme schemes
+        relaxedNonOrthoGaussLaplacianScheme
+        (
+            const fvMesh& mesh,
+            const tmp<surfaceInterpolationScheme<GType>>& igs,
+            const tmp<snGradScheme<Type>>& sngs
+        )
+        :
+            laplacianScheme<Type, GType>(mesh, igs, sngs)
+        {}
+
+
+    //- Destructor
+    virtual ~relaxedNonOrthoGaussLaplacianScheme() = default;
+
+
+    // Member Functions
+
+        static tmp<fvMatrix<Type>> fvmLaplacianUncorrected
+        (
+            const surfaceScalarField& gammaMagSf,
+            const surfaceScalarField& deltaCoeffs,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<fvMatrix<Type>> fvmLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+};
+
+
+// Use macros to emulate partial-specialisation of the Laplacian functions
+// for scalar diffusivity gamma
+
+#define defineFvmLaplacianScalarGamma(Type)                                    \
+                                                                               \
+template<>                                                                     \
+tmp<fvMatrix<Type>>                                                            \
+relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvmLaplacian                \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);                                                                             \
+                                                                               \
+template<>                                                                     \
+tmp<GeometricField<Type, fvPatchField, volMesh>>                               \
+relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvcLaplacian                \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);
+
+
+defineFvmLaplacianScalarGamma(scalar);
+defineFvmLaplacianScalarGamma(vector);
+defineFvmLaplacianScalarGamma(sphericalTensor);
+defineFvmLaplacianScalarGamma(symmTensor);
+defineFvmLaplacianScalarGamma(tensor);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "relaxedNonOrthoGaussLaplacianScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "relaxedNonOrthoGaussLaplacianScheme.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvLaplacianScheme(relaxedNonOrthoGaussLaplacianScheme)
+
+#define declareFvmLaplacianScalarGamma(Type)                                   \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
+Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>::      \
+fvmLaplacian                                                                   \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType;            \
+                                                                               \
+    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf              \
+    (                                                                          \
+        gamma*mesh.magSf()                                                     \
+    );                                                                         \
+                                                                               \
+    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected                         \
+    (                                                                          \
+        gammaMagSf,                                                            \
+        this->tsnGradScheme_().deltaCoeffs(vf),                                \
+        vf                                                                     \
+    );                                                                         \
+    fvMatrix<Type>& fvm = tfvm.ref();                                          \
+                                                                               \
+    if (this->tsnGradScheme_().corrected())                                    \
+    {                                                                          \
+        tmp<SType> tCorr(this->tsnGradScheme_().correction(vf));               \
+        const word corrName(tCorr().name());                                   \
+        tmp<SType> tfaceFluxCorrection(gammaMagSf*tCorr);                      \
+                                                                               \
+        tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));       \
+                                                                               \
+        const word oldName(corrName + "_0");                                   \
+        const scalar relax(vf.mesh().equationRelaxationFactor(corrName));      \
+        const objectRegistry& obr = vf.db();                                   \
+        if (obr.foundObject<SType>(oldName))                                   \
+        {                                                                      \
+            SType& oldCorrection = obr.lookupObjectRef<SType>(oldName);        \
+            trelaxedCorrection.ref() *= relax;                                 \
+            trelaxedCorrection.ref() += (1.0-relax)*oldCorrection;             \
+                                                                               \
+            oldCorrection = trelaxedCorrection();                              \
+        }                                                                      \
+        else                                                                   \
+        {                                                                      \
+            SType* s = new SType(oldName, tfaceFluxCorrection);                \
+            s->store();                                                        \
+        }                                                                      \
+                                                                               \
+        tmp<Field<Type>> tcorr                                                 \
+        (                                                                      \
+            mesh.V()                                                           \
+           *fvc::div                                                           \
+            (                                                                  \
+                trelaxedCorrection()                                           \
+            )().primitiveField()                                               \
+        );                                                                     \
+                                                                               \
+        fvm.source() -= tcorr();                                               \
+                                                                               \
+        if (mesh.fluxRequired(vf.name()))                                      \
+        {                                                                      \
+            fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();            \
+        }                                                                      \
+    }                                                                          \
+                                                                               \
+    return tfvm;                                                               \
+}                                                                              \
+                                                                               \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
+Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>::fvcLaplacian         \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian                \
+    (                                                                          \
+        fvc::div(gamma*this->tsnGradScheme_().snGrad(vf)*mesh.magSf())         \
+    );                                                                         \
+                                                                               \
+    tLaplacian.ref().rename                                                    \
+    (                                                                          \
+        "laplacian(" + gamma.name() + ',' + vf.name() + ')'                    \
+    );                                                                         \
+                                                                               \
+    return tLaplacian;                                                         \
+}
+
+
+declareFvmLaplacianScalarGamma(scalar);
+declareFvmLaplacianScalarGamma(vector);
+declareFvmLaplacianScalarGamma(sphericalTensor);
+declareFvmLaplacianScalarGamma(symmTensor);
+declareFvmLaplacianScalarGamma(tensor);
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pipeCyclic/constant/turbulenceProperties

@@ -24,6 +24,8 @@ RAS
     turbulence      on;
 
     printCoeffs     on;
+
+    implicitKEpsilonSource  true;
 }
 
 

tutorials/incompressible/simpleFoam/pipeCyclic/system/fvSolution

@@ -22,14 +22,14 @@ solvers
         solver          GAMG;
         smoother        GaussSeidel;
 
-        tolerance       1e-06;
+        tolerance       1e-09;
         relTol          0.05;
     }
 
     pFinal
     {
         $p;
-        tolerance       1e-06;
+        tolerance       1e-09;
         relTol          0;
     }
 
@@ -37,14 +37,14 @@ solvers
     {
         solver          smoothSolver;
         smoother        symGaussSeidel;
-        tolerance       1e-05;
+        tolerance       1e-015;
         relTol          0.1;
     }
 
     "(U|k|epsilon)Final"
     {
         $U;
-        tolerance       1e-05;
+        tolerance       1e-15;
         relTol          0;
     }
 }
@@ -61,12 +61,14 @@ PIMPLE
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
+    /*
     residualControl
     {
         p               1e-2;
         U               1e-3;
         "(k|epsilon)"   1e-3;
     }
+    */
 }
 
 relaxationFactors
@@ -79,7 +81,7 @@ relaxationFactors
     {
         U               0.7;
         k               0.7;
-        "epsilon.*"     0.7;
+        epsilon         0.7;
     }
 }
 

tutorials/incompressible/simpleFoam/pitzDaily/constant/turbulenceProperties

@@ -26,6 +26,8 @@ RAS
     turbulence      on;
 
     printCoeffs     on;
+
+    implicitKEpsilonSource  true;
 }