MomentumTransportModels::k-omega model: Improved omega bounding using nutMaxCoeff

Omega lower limit bounding is now based on a maximum turbulence viscosity nut
rather than a minimum omega value which improves stability and robustness of
the k-omega models in case of numerical boundedness problems.

The maximum nut value is calculated by multiplying the laminar viscosity by
nutMaxCoeff which defaults to 1e5 but can be set by the user in the
momentumTransport dictionary.

applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C

@@ -195,7 +195,6 @@ int main(int argc, char *argv[])
                 refCast<const incompressible::RASModel>(turbulence());
 
             omega = (1 - mask)*omega + mask*ce0*sqrt(k)/(Cmu*min(y, ybl));
-            bound(omega, rasModel.omegaMin());
 
             // Do not correct BC - G set by the wall-functions is not available
             // omega.correctBoundaryConditions();

src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.C

@@ -492,6 +492,16 @@ kkLOmega::kkLOmega
             1.17
         )
     ),
+    omegaMin_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "omegaMin",
+            coeffDict_,
+            dimless/dimTime,
+            small
+        )
+    ),
     kt_
     (
         IOobject
@@ -536,12 +546,13 @@ kkLOmega::kkLOmega
             runTime_.name(),
             mesh_
         ),
-        kt_*max(omega_, omegaMin_) + D(kl_) + D(kt_)
+        kt_*omega_
     )
 {
     bound(kt_, kMin_);
     bound(kl_, kMin_);
     bound(omega_, omegaMin_);
+    epsilon_ = kt_*omega_ + D(kl_) + D(kt_);
 
     if (type == typeName)
     {
@@ -751,9 +762,11 @@ void kkLOmega::correct()
     solve(ktEqn);
     bound(kt_, kMin_);
 
+
     // Update total fluctuation kinetic energy dissipation rate
     epsilon_ = kt_*omega_ + Dl + Dt;
 
+
     // Re-calculate turbulent viscosity
     nut_ = nuts + nutl;
     nut_.correctBoundaryConditions();

src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.H

@@ -192,6 +192,9 @@ protected:
             dimensionedScalar Sigmak_;
             dimensionedScalar Sigmaw_;
 
+        //- Lower limit for omega
+        dimensionedScalar omegaMin_;
+
 
         // Fields
 

src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.C

@@ -36,9 +36,17 @@ namespace Foam
 
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
+template<class MomentumTransportModel, class BasicMomentumTransportModel>
+void
+kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::boundOmega()
+{
+    omega_ = max(omega_, k_/(this->nutMaxCoeff_*this->nu()));
+}
+
+
 template<class MomentumTransportModel, class BasicMomentumTransportModel>
 tmp<volScalarField>
-kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::kOmegaSST::F1
+kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::F1
 (
     const volScalarField& CDkOmega
 ) const
@@ -68,8 +76,7 @@ kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::kOmegaSST::F1
 
 template<class MomentumTransportModel, class BasicMomentumTransportModel>
 tmp<volScalarField>
-kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::kOmegaSST::
+kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::F2() const
-F2() const
 {
     tmp<volScalarField> arg2 = min
     (
@@ -86,8 +93,7 @@ F2() const
 
 template<class MomentumTransportModel, class BasicMomentumTransportModel>
 tmp<volScalarField>
-kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::kOmegaSST::
+kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::F3() const
-F3() const
 {
     tmp<volScalarField> arg3 = min
     (
@@ -100,8 +106,7 @@ F3() const
 
 template<class MomentumTransportModel, class BasicMomentumTransportModel>
 tmp<volScalarField>
-kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::kOmegaSST::
+kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::F23() const
-F23() const
 {
     tmp<volScalarField> f23(F2());
 
@@ -127,8 +132,6 @@ void kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::correctNut
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 template<class MomentumTransportModel, class BasicMomentumTransportModel>
 void kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::
 correctNut()
@@ -380,7 +383,7 @@ kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::kOmegaSST
     )
 {
     bound(k_, this->kMin_);
-    bound(omega_, this->omegaMin_);
+    boundOmega();
 }
 
 
@@ -493,7 +496,7 @@ void kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::correct()
         omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
         solve(omegaEqn);
         fvConstraints.constrain(omega_);
-        bound(omega_, this->omegaMin_);
+        boundOmega();
     }
 
     // Turbulent kinetic energy equation
@@ -515,6 +518,7 @@ void kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::correct()
     solve(kEqn);
     fvConstraints.constrain(k_);
     bound(k_, this->kMin_);
+    boundOmega();
 
     correctNut(S2, F23);
 }

src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.H

@@ -200,6 +200,9 @@ protected:
             return blend(F1, gamma1_, gamma2_);
         }
 
+        //- Bound omega
+        void boundOmega();
+
         virtual void correctNut
         (
             const volScalarField& S2,

src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.C

@@ -79,14 +79,14 @@ Foam::LESModel<BasicMomentumTransportModel>::LESModel
         )
     ),
 
-    omegaMin_
+    nutMaxCoeff_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "omegaMin",
+            "nutMaxCoeff",
             LESDict_,
-            dimless/dimTime,
+            dimless,
-            small
+            1e5
         )
     ),
 
@@ -192,6 +192,7 @@ bool Foam::LESModel<BasicMomentumTransportModel>::read()
         delta_().read(LESDict_);
 
         kMin_.readIfPresent(LESDict_);
+        nutMaxCoeff_.readIfPresent(LESDict_);
 
         return true;
     }

src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.H

@@ -82,8 +82,8 @@ protected:
         //- Lower limit of k
         dimensionedScalar kMin_;
 
-        //- Lower limit for omega
+        //- Upper limit coefficient for nut
-        dimensionedScalar omegaMin_;
+        dimensionedScalar nutMaxCoeff_;
 
         //- Run-time selectable generalised Newtonian viscosity model
         autoPtr<laminarModels::generalisedNewtonianViscosityModel>
@@ -176,18 +176,6 @@ public:
             return coeffDict_;
         }
 
-        //- Return the lower allowable limit for k (default: small)
-        const dimensionedScalar& kMin() const
-        {
-            return kMin_;
-        }
-
-        //- Allow kMin to be changed
-        dimensionedScalar& kMin()
-        {
-            return kMin_;
-        }
-
         //- Access function to filter width
         inline const volScalarField& delta() const
         {

src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.C

@@ -79,17 +79,6 @@ Foam::RASModel<BasicMomentumTransportModel>::RASModel
         )
     ),
 
-    omegaMin_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "omegaMin",
-            RASDict_,
-            dimless/dimTime,
-            small
-        )
-    ),
-
     nutMaxCoeff_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -191,7 +180,7 @@ bool Foam::RASModel<BasicMomentumTransportModel>::read()
         coeffDict_ <<= RASDict_.optionalSubDict(type() + "Coeffs");
 
         kMin_.readIfPresent(RASDict_);
-        omegaMin_.readIfPresent(RASDict_);
+        nutMaxCoeff_.readIfPresent(RASDict_);
 
         return true;
     }

src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.H

@@ -75,9 +75,6 @@ protected:
         //- Lower limit of k
         dimensionedScalar kMin_;
 
-        //- Lower limit for omega
-        dimensionedScalar omegaMin_;
-
         //- Upper limit coefficient for nut
         dimensionedScalar nutMaxCoeff_;
 
@@ -163,30 +160,6 @@ public:
         //- Read model coefficients if they have changed
         virtual bool read();
 
-        //- Return the lower allowable limit for k (default: small)
-        const dimensionedScalar& kMin() const
-        {
-            return kMin_;
-        }
-
-        //- Return the lower allowable limit for omega (default: small)
-        const dimensionedScalar& omegaMin() const
-        {
-            return omegaMin_;
-        }
-
-        //- Allow kMin to be changed
-        dimensionedScalar& kMin()
-        {
-            return kMin_;
-        }
-
-        //- Allow omegaMin to be changed
-        dimensionedScalar& omegaMin()
-        {
-            return omegaMin_;
-        }
-
         //- Const access to the coefficients dictionary
         virtual const dictionary& coeffDict() const
         {

src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.C

@@ -37,6 +37,13 @@ namespace RASModels
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class BasicMomentumTransportModel>
+void kOmega<BasicMomentumTransportModel>::boundOmega()
+{
+    omega_ = max(omega_, k_/(this->nutMaxCoeff_*this->nu()));
+}
+
+
 template<class BasicMomentumTransportModel>
 void kOmega<BasicMomentumTransportModel>::correctNut()
 {
@@ -172,7 +179,7 @@ kOmega<BasicMomentumTransportModel>::kOmega
     )
 {
     bound(k_, this->kMin_);
-    bound(omega_, this->omegaMin_);
+    boundOmega();
 
     if (type == typeName)
     {
@@ -260,7 +267,7 @@ void kOmega<BasicMomentumTransportModel>::correct()
     omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
     solve(omegaEqn);
     fvConstraints.constrain(omega_);
-    bound(omega_, this->omegaMin_);
+    boundOmega();
 
 
     // Turbulent kinetic energy equation
@@ -282,6 +289,7 @@ void kOmega<BasicMomentumTransportModel>::correct()
     solve(kEqn);
     fvConstraints.constrain(k_);
     bound(k_, this->kMin_);
+    boundOmega();
 
     correctNut();
 }

src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.H

@@ -96,8 +96,16 @@ protected:
 
     // Protected Member Functions
 
+        //- Bound omega
+        void boundOmega();
+
+        //- Correct the eddy-viscosity nut
         virtual void correctNut();
+
+        //- Source term for the k equation
         virtual tmp<fvScalarMatrix> kSource() const;
+
+        //- Source term for the omega equation
         virtual tmp<fvScalarMatrix> omegaSource() const;
 
 

src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.C

@@ -35,6 +35,15 @@ namespace Foam
 namespace RASModels
 {
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicMomentumTransportModel>
+void kOmega2006<BasicMomentumTransportModel>::boundOmega()
+{
+    omega_ = max(omega_, k_/(this->nutMaxCoeff_*this->nu()));
+}
+
+
 template<class BasicMomentumTransportModel>
 void kOmega2006<BasicMomentumTransportModel>::correctNut
 (
@@ -238,7 +247,7 @@ kOmega2006<BasicMomentumTransportModel>::kOmega2006
     )
 {
     bound(k_, this->kMin_);
-    bound(omega_, this->omegaMin_);
+    boundOmega();
 
     if (type == typeName)
     {
@@ -328,7 +337,7 @@ void kOmega2006<BasicMomentumTransportModel>::correct()
     omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
     solve(omegaEqn);
     fvConstraints.constrain(omega_);
-    bound(omega_, this->omegaMin_);
+    boundOmega();
 
 
     // Turbulent kinetic energy equation
@@ -350,6 +359,7 @@ void kOmega2006<BasicMomentumTransportModel>::correct()
     solve(kEqn);
     fvConstraints.constrain(k_);
     bound(k_, this->kMin_);
+    boundOmega();
 
     correctNut(gradU);
 }

src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.H

@@ -111,8 +111,16 @@ protected:
 
     // Protected Member Functions
 
+        //- Bound omega
+        void boundOmega();
+
+        //- Correct the eddy-viscosity nut
         virtual void correctNut();
+
+        //- Source term for the k equation
         virtual tmp<fvScalarMatrix> kSource() const;
+
+        //- Source term for the omega equation
         virtual tmp<fvScalarMatrix> omegaSource() const;