From d3d00555adcc8d9dd52c6de806c9c52753ef9b54 Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Fri, 22 Jun 2018 10:50:02 +0200
Subject: [PATCH] [freeflow][komega] Add dissipation limiter
---
dumux/freeflow/rans/twoeq/komega/model.hh | 3 ---
dumux/freeflow/rans/twoeq/komega/problem.hh | 2 +-
.../twoeq/komega/staggered/localresidual.hh | 8 ++++---
.../rans/twoeq/komega/volumevariables.hh | 21 ++++++++++++-------
4 files changed, 20 insertions(+), 14 deletions(-)
diff --git a/dumux/freeflow/rans/twoeq/komega/model.hh b/dumux/freeflow/rans/twoeq/komega/model.hh
index 52ab38eaca..eb300623f1 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 f7b1651447..d5a4c6704b 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 5d86d36043..8d30b51e06 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 aad519d26c..98626dcc72 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();
}
/*!
--
GitLab