diff --git a/doc/handbook/dumux-handbook.bib b/doc/handbook/dumux-handbook.bib
index 4bfc36d0e30c0cceb43b464df2b4b2d909ff3b25..6e14fbf946e5d1a411e815d1c485964f6865787a 100644
--- a/doc/handbook/dumux-handbook.bib
+++ b/doc/handbook/dumux-handbook.bib
@@ -1714,3 +1714,14 @@ url={http://dx.doi.org/10.1007/s11242-015-0599-1}
doi = {10.2514/1.36541},
url = {https://doi.org/10.2514/1.36541}
}
+
+@Book{Versteeg2009a,
+ title = {{An Introduction to Computational Fluid Dynamics}},
+ publisher = {Pearson Education},
+ year = {2009},
+ author = {Versteeg, Henk and Malalasekra, Weeratunge},
+ isbn = {978-0-13-127498-3},
+ pages = {XII, 503},
+ address = {Harlow},
+ edition = {2},
+}
diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh
index c1e9d532414cfee8a603a20fd1f224261d575b50..18f14136f434a88cdfd71b4a8aadc6adfe0d213c 100644
--- a/dumux/freeflow/navierstokes/problem.hh
+++ b/dumux/freeflow/navierstokes/problem.hh
@@ -76,6 +76,7 @@ class NavierStokesProblem : public NavierStokesParentProblem<TypeTag>
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
enum {
@@ -190,6 +191,14 @@ public:
const SubControlVolumeFace& localSubFace) const
{ return FacePrimaryVariables(0.0); }
+ //! \brief Returns an additional wall function flux for cell-centered quantities (only needed for RANS models)
+ CellCenterPrimaryVariables wallFunction(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ { return CellCenterPrimaryVariables(0.0); }
+
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 7042b56ff4ac25ba019d9f7312218ef203781752..2a96c338f16d62dd297923b5463f529459faaa4a 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -79,6 +79,9 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();
+ static constexpr bool isCompositional = ModelTraits::numComponents() > 1;
+
public:
using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance(), (ModelTraits::numComponents() > 1)>;
@@ -222,6 +225,7 @@ protected:
if (scvf.boundary())
{
const auto bcTypes = problem.boundaryTypes(element, scvf);
+ const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
// treat Dirichlet and outflow BCs
FluxVariables fluxVars;
@@ -231,19 +235,27 @@ protected:
EnergyLocalResidual::heatFlux(boundaryFlux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
// treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
+ static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
if(bcTypes.hasNeumann())
{
- static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
{
if(bcTypes.isNeumann(eqIdx + cellCenterOffset))
{
- const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx + cellCenterOffset]
- * extrusionFactor * scvf.area();
+ * extrusionFactor * scvf.area();
}
}
}
+ for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+ {
+ // use a wall function
+ if(problem.useWallFunction(element, scvf, eqIdx + cellCenterOffset))
+ {
+ boundaryFlux[eqIdx] = problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx]
+ * extrusionFactor * scvf.area();
+ }
+ }
residual += boundaryFlux;
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 5b703a38888460849c5d9a31ed7b25775080ca77..3f7c59cf577c8ee40a684308d91bb61632416dfc 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -68,9 +68,20 @@ class KEpsilonProblem : public RANSProblem<TypeTag>
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
- using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ using Indices = typename ModelTraits::Indices;
+
+ static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();
+ static constexpr bool isCompositional = ModelTraits::numComponents() > 1;
+
+ // account for the offset of the cell center privars within the PrimaryVariables container
+ static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;
+ static_assert(cellCenterOffset == ModelTraits::dim(), "cellCenterOffset must equal dim for staggered NavierStokes");
public:
+ static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
//! The constructor sets the gravity, if desired by the user.
KEpsilonProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
@@ -319,6 +330,119 @@ public:
* elemVolVars[scvf.insideScvIdx()].density());
}
+ //! \brief Returns the flux for non-isothermal and compositional RANS models
+ template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+ typename std::enable_if_t<eB && compositional, int> = 0>
+ CellCenterPrimaryVariables wallFunction(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ return wallFunctionComponent(element, fvGeometry, elemVolVars, elemFaceVars, scvf)
+ + wallFunctionEnergy(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+ }
+
+ //! \brief Returns the flux for isothermal and compositional RANS models
+ template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+ typename std::enable_if_t<!eB && compositional, int> = 0>
+ CellCenterPrimaryVariables wallFunction(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ { return wallFunctionComponent(element, fvGeometry, elemVolVars, elemFaceVars, scvf); }
+
+ //! \brief Returns the flux for non-isothermal RANS models
+ template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+ typename std::enable_if_t<eB && !compositional, int> = 0>
+ CellCenterPrimaryVariables wallFunction(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ { return wallFunctionEnergy(element, fvGeometry, elemVolVars, elemFaceVars, scvf); }
+
+ //! \brief Returns the flux for isothermal RANS models
+ template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+ typename std::enable_if_t<!eB && !compositional, int> = 0>
+ CellCenterPrimaryVariables wallFunction(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ { return CellCenterPrimaryVariables(0.0); }
+
+ //! \brief Returns the component wall-function flux
+ CellCenterPrimaryVariables wallFunctionComponent(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ auto wallFunctionFlux = CellCenterPrimaryVariables(0.0);
+ unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element);
+
+ // component mass fluxes
+ for (int compIdx = 0; compIdx < ModelTraits::numComponents(); ++compIdx)
+ {
+ if (Indices::replaceCompEqIdx == compIdx)
+ continue;
+
+ Scalar schmidtNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity()
+ / elemVolVars[scvf.insideScvIdx()].diffusionCoefficient(compIdx);
+ Scalar massConversionFactor = useMoles ? 1.0
+ : FluidSystem::molarMass(compIdx);
+ wallFunctionFlux[compIdx] +=
+ -1.0 * (asImp_().dirichlet(element, scvf)[Indices::conti0EqIdx + compIdx]
+ - elemVolVars[scvf.insideScvIdx()].moleFraction(compIdx))
+ * elemVolVars[scvf.insideScvIdx()].molarDensity()
+ * uStarNominal(elementID)
+ / asImp_().turbulentSchmidtNumber()
+ / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementID) * 9.793)
+ + pFunction(schmidtNumber, asImp_().turbulentSchmidtNumber()));
+ }
+
+ return wallFunctionFlux;
+ }
+
+ //! \brief Returns the energy wall-function flux
+ CellCenterPrimaryVariables wallFunctionEnergy(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ auto wallFunctionFlux = CellCenterPrimaryVariables(0.0);
+ unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element);
+ // energy fluxes
+ Scalar prandtlNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity()
+ * elemVolVars[scvf.insideScvIdx()].density()
+ * elemVolVars[scvf.insideScvIdx()].heatCapacity()
+ / elemVolVars[scvf.insideScvIdx()].thermalConductivity();
+ wallFunctionFlux[Indices::energyBalanceIdx - cellCenterOffset] +=
+ -1.0 * (asImp_().dirichlet(element, scvf)[Indices::temperatureIdx]
+ - elemVolVars[scvf.insideScvIdx()].temperature())
+ * elemVolVars[scvf.insideScvIdx()].density()
+ * elemVolVars[scvf.insideScvIdx()].heatCapacity()
+ * uStarNominal(elementID)
+ / asImp_().turbulentPrandtlNumber()
+ / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementID) * 9.793)
+ + pFunction(prandtlNumber, asImp_().turbulentPrandtlNumber()));
+
+ return wallFunctionFlux;
+ }
+
+ //! \brief Returns the value of the P-function after Jayatilleke \cite Versteeg2009a
+ const Scalar pFunction(Scalar molecularNumber, Scalar turbulentNumber) const
+ {
+ using std::pow;
+ using std::exp;
+ return 9.24
+ * (pow(molecularNumber / turbulentNumber, 0.75) - 1.0)
+ * (1.0 + 0.28 * exp(-0.007 * molecularNumber / turbulentNumber));
+ }
+
//! \brief Returns the \$f C_{\mu} \$f constant
const Scalar cMu() const
{