src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega: Corrected errors in implementation and from original paper according to

Furst, J. (2013).
    Numerical simulation of transitional flows with laminar kinetic energy.
    Engineering MECHANICS, 20(5), 379-388.

Thanks to Jan-Niklas Klatt for analysing problems with and correcting
the implementation and testing corrections to the model proposed by
Furst.

src/TurbulenceModels/incompressible/incompressibleTurbulenceModel.H

@@ -61,6 +61,9 @@ protected:
 
         geometricOneField rho_;
 
+
+    // Protected member functions
+
         //- ***HGW Temporary function to be removed when the run-time selectable
         //  thermal transport layer is complete
         virtual void correctNut()

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.C

@@ -56,16 +56,16 @@ tmp<volScalarField> kkLOmega::fINT() const
     (
         min
         (
-            kl_/(Cint_*(kl_ + kt_ + kMin_)),
+            kt_/(Cint_*(kl_ + kt_ + kMin_)),
             dimensionedScalar("1.0", dimless, 1.0)
         )
     );
 }
 
 
-tmp<volScalarField> kkLOmega::fSS(const volScalarField& omega) const
+tmp<volScalarField> kkLOmega::fSS(const volScalarField& Omega) const
 {
-    return(exp(-sqr(Css_*nu()*omega/(kt_ + kMin_))));
+    return(exp(-sqr(Css_*nu()*Omega/(kt_ + kMin_))));
 }
 
 
@@ -75,16 +75,17 @@ tmp<volScalarField> kkLOmega::Cmu(const volScalarField& S) const
 }
 
 
-tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& Rew) const
+tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& ReOmega) const
 {
-    return(scalar(1) - exp(-sqr(max(Rew - CtsCrit_, scalar(0)))/Ats_));
+    return(scalar(1) - exp(-sqr(max(ReOmega - CtsCrit_, scalar(0)))/Ats_));
 }
 
 
 tmp<volScalarField> kkLOmega::fTaul
 (
     const volScalarField& lambdaEff,
-    const volScalarField& ktL
+    const volScalarField& ktL,
+    const volScalarField& Omega
 ) const
 {
     return
@@ -97,7 +98,7 @@ tmp<volScalarField> kkLOmega::fTaul
             (
                 sqr
                 (
-                    lambdaEff*omega_
+                    lambdaEff*Omega
                   + dimensionedScalar
                     (
                         "ROOTVSMALL",
@@ -133,8 +134,8 @@ tmp<volScalarField> kkLOmega::fOmega
         scalar(1)
       - exp
         (
-            -0.41
+           -0.41
-           * pow4
+           *pow4
             (
                 lambdaEff
               / (
@@ -152,7 +153,7 @@ tmp<volScalarField> kkLOmega::fOmega
 }
 
 
-tmp<volScalarField> kkLOmega::gammaBP(const volScalarField& omega) const
+tmp<volScalarField> kkLOmega::phiBP(const volScalarField& Omega) const
 {
     return
     (
@@ -162,11 +163,11 @@ tmp<volScalarField> kkLOmega::gammaBP(const volScalarField& omega) const
             (
                 kt_/nu()
              / (
-                    omega
+                    Omega
                   + dimensionedScalar
                     (
                         "ROTVSMALL",
-                        omega.dimensions(),
+                        Omega.dimensions(),
                         ROOTVSMALL
                     )
                 )
@@ -179,7 +180,7 @@ tmp<volScalarField> kkLOmega::gammaBP(const volScalarField& omega) const
 }
 
 
-tmp<volScalarField> kkLOmega::gammaNAT
+tmp<volScalarField> kkLOmega::phiNAT
 (
     const volScalarField& ReOmega,
     const volScalarField& fNatCrit
@@ -200,12 +201,19 @@ tmp<volScalarField> kkLOmega::gammaNAT
 }
 
 
+tmp<volScalarField> kkLOmega::D(const volScalarField& k) const
+{
+    return nu()*magSqr(fvc::grad(sqrt(k)));
+}
+
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void kkLOmega::correctNut()
 {
-    nut_ = kt_/omega_;
+    // Currently this function is not implemented due to the complexity of
-    nut_.correctBoundaryConditions();
+    // evaluating nut.  Better calculate nut at the end of correct()
+    notImplemented("kkLOmega::correctNut()");
 }
 
 
@@ -490,18 +498,6 @@ kkLOmega::kkLOmega
         ),
         mesh_
     ),
-    omega_
-    (
-        IOobject
-        (
-            IOobject::groupName("omega", U.group()),
-            runTime_.timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
     kl_
     (
         IOobject
@@ -514,17 +510,38 @@ kkLOmega::kkLOmega
         ),
         mesh_
     ),
-    y_(wallDist::New(mesh_).y())
+    omega_
+    (
+        IOobject
+        (
+            IOobject::groupName("omega", U.group()),
+            runTime_.timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    epsilon_
+    (
+        IOobject
+        (
+            "epsilon",
+            runTime_.timeName(),
+            mesh_
+        ),
+        kt_*omega_ + D(kl_) + D(kt_)
+    ),    y_(wallDist::New(mesh_).y())
 {
     bound(kt_, kMin_);
     bound(kl_, kMin_);
     bound(omega_, omegaMin_);
+    bound(epsilon_, epsilonMin_);
 
     if (type == typeName)
     {
-        // This is only an approximate nut, so only good for initialization
+        // Evaluating nut_ is complex so start from the field read from file
-        // not restart
+        nut_.correctBoundaryConditions();
-        // correctNut();
 
         printCoeffs(type);
     }
@@ -584,22 +601,28 @@ void kkLOmega::correct()
         return;
     }
 
+    const volScalarField lambdaT(sqrt(kt_)/(omega_ + omegaMin_));
 
-    const volScalarField kT(kt_ + kl_);
-
-    const volScalarField lambdaT(sqrt(kT)/(omega_ + omegaMin_));
     const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
 
     const volScalarField fw
     (
-        lambdaEff/(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL))
+        pow
+        (
+            lambdaEff
+           /(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL)),
+            2.0/3.0
+        )
     );
 
-    const volTensorField gradU(fvc::grad(U_));
+    tmp<volTensorField> tgradU(fvc::grad(U_));
-    const volScalarField omega(sqrt(2.0)*mag(skew(gradU)));
+    const volTensorField& gradU = tgradU();
-    const volScalarField S2(2.0*magSqr(symm(gradU)));
 
-    const volScalarField ktS(fSS(omega)*fw*kt_);
+    const volScalarField Omega(sqrt(2.0)*mag(skew(gradU)));
+
+    const volScalarField S2(2.0*magSqr(dev(symm(gradU))));
+
+    const volScalarField ktS(fSS(Omega)*fw*kt_);
 
     const volScalarField nuts
     (
@@ -607,20 +630,19 @@ void kkLOmega::correct()
        *fINT()
        *Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
     );
-
     const volScalarField Pkt(nuts*S2);
 
     const volScalarField ktL(kt_ - ktS);
-    const volScalarField ReOmega(sqr(y_)*omega/nu());
+    const volScalarField ReOmega(sqr(y_)*Omega/nu());
     const volScalarField nutl
     (
         min
         (
-            C11_*fTaul(lambdaEff, ktL)*omega*sqr(lambdaEff)
+            C11_*fTaul(lambdaEff, ktL, Omega)*Omega*sqr(lambdaEff)
-          * sqrt(ktL)*lambdaEff/nu()
+           *sqrt(ktL)*lambdaEff/nu()
-          + C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*omega
+          + C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*Omega
-          ,
+        ,
-            0.5*(kl_ + ktL)/sqrt(S2)
+            0.5*(kl_ + ktL)/(sqrt(S2) + omegaMin_)
         )
     );
 
@@ -637,60 +659,18 @@ void kkLOmega::correct()
 
     const volScalarField Rbp
     (
-        CR_*(1.0 - exp(-gammaBP(omega)()/Abp_))*omega_
+        CR_*(1.0 - exp(-phiBP(Omega)()/Abp_))*omega_
-      / (fw + fwMin)
+       /(fw + fwMin)
     );
 
     const volScalarField fNatCrit(1.0 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
+
     // Natural source term divided by kl_
     const volScalarField Rnat
     (
-        CrNat_*(1.0 - exp(-gammaNAT(ReOmega, fNatCrit)/Anat_))*omega
+        CrNat_*(1.0 - exp(-phiNAT(ReOmega, fNatCrit)/Anat_))*Omega
     );
 
-    const volScalarField Dt(nu()*magSqr(fvc::grad(sqrt(kt_))));
-
-    // Turbulent kinetic energy equation
-    tmp<fvScalarMatrix> ktEqn
-    (
-        fvm::ddt(kt_)
-      + fvm::div(phi_, kt_)
-      - fvm::laplacian(DkEff(alphaTEff), kt_, "laplacian(alphaTEff,kt)")
-     ==
-        Pkt
-      + (Rbp + Rnat)*kl_
-      - Dt
-      - fvm::Sp(omega_, kt_)
-    );
-
-    ktEqn().relax();
-    ktEqn().boundaryManipulate(kt_.boundaryField());
-
-    solve(ktEqn);
-    bound(kt_, kMin_);
-
-
-    const volScalarField Dl(nu()*magSqr(fvc::grad(sqrt(kl_))));
-
-    // Laminar kinetic energy equation
-    tmp<fvScalarMatrix> klEqn
-    (
-        fvm::ddt(kl_)
-      + fvm::div(phi_, kl_)
-      - fvm::laplacian(nu(), kl_, "laplacian(nu,kl)")
-     ==
-        Pkl
-      - fvm::Sp(Rbp, kl_)
-      - fvm::Sp(Rnat, kl_)
-      - Dl
-    );
-
-    klEqn().relax();
-    klEqn().boundaryManipulate(kl_.boundaryField());
-
-    solve(klEqn);
-    bound(kl_, kMin_);
-
 
     omega_.boundaryField().updateCoeffs();
 
@@ -699,33 +679,74 @@ void kkLOmega::correct()
     (
         fvm::ddt(omega_)
       + fvm::div(phi_, omega_)
-      - fvm::laplacian
+      - fvm::laplacian(DomegaEff(alphaTEff), omega_)
-        (
-            DomegaEff(alphaTEff),
-            omega_,
-            "laplacian(alphaTEff,omega)"
-        )
      ==
         Cw1_*Pkt*omega_/(kt_ + kMin_)
-      + fvm::SuSp
+      - fvm::SuSp
         (
-            (CwR_/(fw + fwMin) - 1.0)*kl_*(Rbp + Rnat)/(kt_ + kMin_)
+            (1.0 - CwR_/(fw + fwMin))*kl_*(Rbp + Rnat)/(kt_ + kMin_)
           , omega_
         )
-      - fvm::Sp(Cw2_*omega_, omega_)
+      - fvm::Sp(Cw2_*sqr(fw)*omega_, omega_)
       + (
             Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)
         )().dimensionedInternalField()/pow3(y_.dimensionedInternalField())
     );
 
-
     omegaEqn().relax();
     omegaEqn().boundaryManipulate(omega_.boundaryField());
 
     solve(omegaEqn);
     bound(omega_, omegaMin_);
 
-    // Re-calculate viscosity
+
+    const volScalarField Dl(D(kl_));
+
+    // Laminar kinetic energy equation
+    tmp<fvScalarMatrix> klEqn
+    (
+        fvm::ddt(kl_)
+      + fvm::div(phi_, kl_)
+      - fvm::laplacian(nu(), kl_)
+     ==
+        Pkl
+      - fvm::Sp(Rbp + Rnat + Dl/(kl_ + kMin_), kl_)
+    );
+
+    klEqn().relax();
+    klEqn().boundaryManipulate(kl_.boundaryField());
+
+    solve(klEqn);
+    bound(kl_, kMin_);
+
+
+    const volScalarField Dt(D(kt_));
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> ktEqn
+    (
+        fvm::ddt(kt_)
+      + fvm::div(phi_, kt_)
+      - fvm::laplacian(DkEff(alphaTEff), kt_)
+     ==
+        Pkt
+      + (Rbp + Rnat)*kl_
+      - fvm::Sp(omega_ + Dt/(kt_+ kMin_), kt_)
+    );
+
+    ktEqn().relax();
+    ktEqn().boundaryManipulate(kt_.boundaryField());
+
+    solve(ktEqn);
+    bound(kt_, kMin_);
+
+
+    // Update total fluctuation kinetic energy dissipation rate
+    epsilon_ = kt_*omega_ + Dl + Dt;
+    bound(epsilon_, epsilonMin_);
+
+
+    // Re-calculate turbulent viscosity
     nut_ = nuts + nutl;
     nut_.correctBoundaryConditions();
 }

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/kkLOmega/kkLOmega.H

@@ -31,7 +31,7 @@ Description
     Low Reynolds-number k-kl-omega turbulence model for
     incompressible flows.
 
-    Turbulence model described in:
+    This turbulence model is described in:
     \verbatim
         Walters, D. K., & Cokljat, D. (2008).
         A three-equation eddy-viscosity model for Reynolds-averaged
@@ -39,6 +39,17 @@ Description
         Journal of Fluids Engineering, 130(12), 121401.
     \endverbatim
 
+    however the paper contains several errors which must be corrected for the
+    model to operation correctly as explained in
+
+    \verbatim
+        Furst, J. (2013).
+        Numerical simulation of transitional flows with laminar kinetic energy.
+        Engineering MECHANICS, 20(5), 379-388.
+    \endverbatim
+
+    All these corrections and updates are included in this implementation.
+
     The default model coefficients are
     \verbatim
         kkLOmegaCoeffs
@@ -116,7 +127,8 @@ class kkLOmega
         tmp<volScalarField> fTaul
         (
             const volScalarField& lambdaEff,
-            const volScalarField& ktL
+            const volScalarField& ktL,
+            const volScalarField& omega
         ) const;
 
         tmp<volScalarField> alphaT
@@ -132,14 +144,16 @@ class kkLOmega
             const volScalarField& lambdaT
         ) const;
 
-        tmp<volScalarField> gammaBP(const volScalarField& omega) const;
+        tmp<volScalarField> phiBP(const volScalarField& omega) const;
 
-        tmp<volScalarField> gammaNAT
+        tmp<volScalarField> phiNAT
         (
             const volScalarField& ReOmega,
             const volScalarField& fNatCrit
         ) const;
 
+        tmp<volScalarField> D(const volScalarField& k) const;
+
 
 protected:
 
@@ -179,8 +193,9 @@ protected:
         // Fields
 
             volScalarField kt_;
-            volScalarField omega_;
             volScalarField kl_;
+            volScalarField omega_;
+            volScalarField epsilon_;
 
             //- Wall distance
             //  Note: different to wall distance in parent RASModel
@@ -250,7 +265,7 @@ public:
         }
 
         //- Return the turbulence kinetic energy
-        virtual tmp<volScalarField> k() const
+        virtual tmp<volScalarField> kt() const
         {
             return kt_;
         }
@@ -261,8 +276,8 @@ public:
             return omega_;
         }
 
-        //- Return the turbulence kinetic energy dissipation rate
+        //- Return the total fluctuation kinetic energy
-        virtual tmp<volScalarField> epsilon() const
+        virtual tmp<volScalarField> k() const
         {
             return tmp<volScalarField>
             (
@@ -270,16 +285,22 @@ public:
                 (
                     IOobject
                     (
-                        "epsilon",
+                        "k",
                         mesh_.time().timeName(),
                         mesh_
                     ),
-                    kt_*omega_,
+                    kt_ + kl_,
                     omega_.boundaryField().types()
                 )
             );
         }
 
+        //- Return the total fluctuation kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return epsilon_;
+        }
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 };