[freeflow][komega] Add dissipation limiter

dumux/freeflow/rans/twoeq/komega/model.hh

@@ -104,9 +104,6 @@ struct KOmegaModelTraits : RANSModelTraits<dimension>
 
     //! The indices
     using Indices = KOmegaIndices<dim(), numComponents()>;
-
-    //! The model includes a limiter to the production term
-    static constexpr bool enableKOmegaProductionLimiter() { return true; }
 };
 
 ///////////////////////////////////////////////////////////////////////////

dumux/freeflow/rans/twoeq/komega/problem.hh

@@ -115,7 +115,7 @@ public:
                 // NOTE: then update the volVars
                 VolumeVariables volVars;
                 volVars.update(elemSol, asImp_(), element, scv);
-                storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity();
+                storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity(*this);
             }
         }
     }

dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh

@@ -103,11 +103,13 @@ public:
         const auto& volVars = elemVolVars[scv];
 
         // production
+        static const auto enableKOmegaProductionLimiter
+            = getParamFromGroup<bool>(problem.paramGroup(), "KOmega.EnableProductionLimiter", false);
         Scalar productionTerm = 2.0 * volVars.kinematicEddyViscosity() * volVars.stressTensorScalarProduct();
-        if (ModelTraits::enableKOmegaProductionLimiter())
+        if (enableKOmegaProductionLimiter)
         {
-          Scalar productionAlternative = 20.0 * volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
-          productionTerm = min(productionTerm, productionAlternative);
+            Scalar productionAlternative = 20.0 * volVars.betaK() * volVars.turbulentKineticEnergy() * volVars.dissipation();
+            productionTerm = min(productionTerm, productionAlternative);
         }
         source[turbulentKineticEnergyEqIdx] += productionTerm;
         source[dissipationEqIdx] += volVars.alpha() * volVars.dissipation() / volVars.turbulentKineticEnergy() * productionTerm;

dumux/freeflow/rans/twoeq/komega/volumevariables.hh

@@ -47,6 +47,7 @@ class KOmegaVolumeVariables
     using NavierStokesParentType = NSVolumeVariables;
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
+    using ModelTraits = typename Traits::ModelTraits;
 
     static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
     static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
@@ -101,23 +102,29 @@ public:
         if (problem.useStoredEddyViscosity_)
             dynamicEddyViscosity_ = problem.storedDynamicEddyViscosity_[RANSParentType::elementID()];
         else
-            dynamicEddyViscosity_ = calculateEddyViscosity();
+            dynamicEddyViscosity_ = calculateEddyViscosity(problem);
         calculateEddyDiffusivity(problem);
     }
 
     /*!
      * \brief Returns the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$.
      */
-    Scalar calculateEddyViscosity()
+    template<class Problem>
+    Scalar calculateEddyViscosity(const Problem& problem)
     {
         using std::sqrt;
         using std::max;
 
-        //! Use the Dissipation limiter proposed in wilcox08
-        Scalar limitiedDissipation = (7.0 / 8.0) * sqrt( 2.0 * stressTensorScalarProduct()
-                                     * stressTensorScalarProduct() / betaK() );
-        Scalar Wbar = max(dissipation(), limitiedDissipation);
-        return turbulentKineticEnergy() / Wbar * NavierStokesParentType::density();
+        // use the Dissipation limiter proposed in wilcox2008
+        Scalar limitiedDissipation = std::numeric_limits<Scalar>::min();
+        static const auto enableKOmegaDissipationLimiter
+            = getParamFromGroup<bool>(problem.paramGroup(), "KOmega.EnableDissipationLimiter", true);
+        if (enableKOmegaDissipationLimiter)
+        {
+            limitiedDissipation = (7.0 / 8.0) * sqrt(2.0 * stressTensorScalarProduct() / betaK());
+        }
+        return turbulentKineticEnergy() / max(dissipation(), limitiedDissipation)
+               * NavierStokesParentType::density();
     }
 
     /*!