diff --git a/dumux/freeflow/rans/twoeq/komega/model.hh b/dumux/freeflow/rans/twoeq/komega/model.hh index 52ab38eaca835869eab8a008f2c47dea9c87508f..eb300623f145b4a474ec94a30e12c128808b0684 100644 --- a/dumux/freeflow/rans/twoeq/komega/model.hh +++ b/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; } }; /////////////////////////////////////////////////////////////////////////// diff --git a/dumux/freeflow/rans/twoeq/komega/problem.hh b/dumux/freeflow/rans/twoeq/komega/problem.hh index f7b1651447284852e14d3e46f23735649b91a596..d5a4c6704b7d0b6dd0f72e8a18140c810d4352c1 100644 --- a/dumux/freeflow/rans/twoeq/komega/problem.hh +++ b/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); } } } diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh index 5d86d36043d4eef25b023506671f8944759d5007..8d30b51e06d8b6308556473a05cbae97178a61b2 100644 --- a/dumux/freeflow/rans/twoeq/komega/staggered/localresidual.hh +++ b/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; diff --git a/dumux/freeflow/rans/twoeq/komega/volumevariables.hh b/dumux/freeflow/rans/twoeq/komega/volumevariables.hh index aad519d26c5f3a12af949552d2653784b46eef26..98626dcc72ac7deb0f14a799c1ba0da0f0b83988 100644 --- a/dumux/freeflow/rans/twoeq/komega/volumevariables.hh +++ b/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(); } /*!