diff --git a/doc/doxygen/modules.txt b/doc/doxygen/modules.txt
index 475280dc7e6e65ccac1e1523a08f27c5e3f3f5aa..3958a6f442b9083c1cd0aef39e4a598c4776e3ab 100644
--- a/doc/doxygen/modules.txt
+++ b/doc/doxygen/modules.txt
@@ -196,6 +196,12 @@
* \brief Two-equation turbulence models
* \copydetails ./freeflow/rans/model.hh
*/
+ /*!
+ * \ingroup TwoEqModel
+ * \defgroup KEpsilonModel K-epsilon model
+ * \brief K-epsilon model
+ * \copydetails ./freeflow/rans/twoeq/kepsilon/model.hh
+ */
/*!
* \ingroup TwoEqModel
* \defgroup LowReKEpsilonModel Low-Re k-epsilon model
diff --git a/dumux/freeflow/rans/twoeq/CMakeLists.txt b/dumux/freeflow/rans/twoeq/CMakeLists.txt
index 445c7e40a012364eff119bcf5fe946d64761f461..ee8d7de6fc7399445c4ec478c02a4a856e8b926d 100644
--- a/dumux/freeflow/rans/twoeq/CMakeLists.txt
+++ b/dumux/freeflow/rans/twoeq/CMakeLists.txt
@@ -1 +1,2 @@
+add_subdirectory("kepsilon")
add_subdirectory("lowrekepsilon")
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/CMakeLists.txt b/dumux/freeflow/rans/twoeq/kepsilon/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1f77d858e5177b169536da8dca85bb9d9f7052cb
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/CMakeLists.txt
@@ -0,0 +1,11 @@
+#install headers
+install(FILES
+fluxvariables.hh
+indices.hh
+localresidual.hh
+model.hh
+# models.hh
+problem.hh
+volumevariables.hh
+vtkoutputfields.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/twoeq/kepsilon)
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/fluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ec1000024cefe7f6ae7190e203423d88153100cf
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/fluxvariables.hh
@@ -0,0 +1,49 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+ /*!
+ * \file
+ * \ingroup KEpsilonModel
+ * \copydoc Dumux::KEpsilonFluxVariables
+ */
+#ifndef DUMUX_KEPSILON_FLUXVARIABLES_HH
+#define DUMUX_KEPSILON_FLUXVARIABLES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template<class TypeTag, class BaseFluxVariables, DiscretizationMethod discMethod>
+class KEpsilonFluxVariablesImpl;
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief The flux variables class for the k-epsilon model.
+ This is a convenience alias for that actual,
+ discretization-specific flux variables.
+ * \note Not all specializations are currently implemented
+ */
+template<class TypeTag, class BaseFluxVariables>
+using KEpsilonFluxVariables = KEpsilonFluxVariablesImpl<TypeTag, BaseFluxVariables, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/indices.hh b/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..98a2fb402e6950dcd15df878a8df6d75dec41180
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
@@ -0,0 +1,52 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup KEpsilonModel
+ * \copydoc Dumux::KEpsilonIndices
+ */
+#ifndef DUMUX_KEPSILON_INDICES_HH
+#define DUMUX_KEPSILON_INDICES_HH
+
+#include <dumux/freeflow/navierstokes/indices.hh>
+
+namespace Dumux {
+
+// \{
+/*!
+ * \ingroup KEpsilonModel
+ * \brief The common indices for the isothermal k-epsilon model.
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam numComponents The number of considered transported components
+ */
+template<int dimension, int numComponents>
+struct KEpsilonIndices : public NavierStokesIndices<dimension>
+{
+public:
+ static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;
+ static constexpr auto turbulentKineticEnergyIdx = turbulentKineticEnergyEqIdx;
+ static constexpr auto dissipationEqIdx = turbulentKineticEnergyEqIdx + 1;
+ static constexpr auto dissipationIdx = dissipationEqIdx;
+};
+
+// \}
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/localresidual.hh b/dumux/freeflow/rans/twoeq/kepsilon/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ac340be38d7e75eebef787d507a23b2ac2ef4e51
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/localresidual.hh
@@ -0,0 +1,51 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+ /*!
+ * \file
+ * \ingroup KEpsilonModel
+ * \copydoc Dumux::KEpsilonResidual
+ */
+#ifndef DUMUX_KEPSILON_LOCAL_RESIDUAL_HH
+#define DUMUX_KEPSILON_LOCAL_RESIDUAL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/localresidual.hh>
+#include <dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template<class TypeTag, class BaseLocalResidual, DiscretizationMethod discMethod>
+class KEpsilonResidualImpl;
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief The local residual class for the k-epsilon model.
+ This is a convenience alias for the actual,
+ discretization-specific local residual.
+ * \note Not all specializations are currently implemented
+ */
+template<class TypeTag, class BaseLocalResidual>
+using KEpsilonResidual = KEpsilonResidualImpl<TypeTag, BaseLocalResidual, GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod>;
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/model.hh b/dumux/freeflow/rans/twoeq/kepsilon/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..94f28bb80e95c47c06da955593dfaa59d3c3ae18
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/model.hh
@@ -0,0 +1,227 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup KEpsilonModel
+ *
+ * \brief A single-phase, isothermal k-epsilon model
+ *
+ * \copydoc RANSModel
+ *
+ * The k-epsilon models calculate the eddy viscosity with two additional PDEs,
+ * one for the turbulent kinetic energy (k) and for the dissipation (\f$ \varepsilon \f$).
+ * \todo fix everything below here
+ * The model uses the one proposed by Chien \cite Chien1982a.
+ * A good overview and additional models are given in Patel et al. \cite Patel1985a.
+ *
+ * The turbulent kinetic energy balance is identical with the one from the k-epsilon model,
+ * but the dissipation includes a dampening function (\f$ D_\varepsilon \f$):
+ * \f$ \varepsilon = \tilde{\varepsilon} + D_\varepsilon \f$:
+ *
+ * \f[
+ * \frac{\partial \left( k \right)}{\partial t}
+ * + \nabla \cdot \left( \textbf{v} k \right)
+ * - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_\text{k}} \right) \nabla k \right)
+ * - 2 \nu_\text{t} \textbf{S} \cdot \textbf{S}
+ * + \tilde{\varepsilon}
+ * + D_\varepsilon
+ * = 0
+ * \f].
+ *
+ * The dissipation balance is changed by introducing additional functions
+ * (\f$ E_\text{k}\f$, \f$ f_1 \f$, and \f$ f_2 \f$) to account for a dampening towards the wall:
+ * \f[
+ * \frac{\partial \left( \tilde{\varepsilon} \right)}{\partial t}
+ * + \nabla \cdot \left( \textbf{v} \tilde{\varepsilon} \right)
+ * - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_{\varepsilon}} \right) \nabla \tilde{\varepsilon} \right)
+ * - C_{1\tilde{\varepsilon}} f_1 \frac{\tilde{\varepsilon}}{k} 2 \nu_\text{t} \textbf{S} \cdot \textbf{S}
+ * + C_{2\tilde{\varepsilon}} f_2 \frac{\tilde{\varepsilon}^2}{k}
+ * - E_\text{k}
+ * = 0
+ * \f].
+ *
+ * The kinematic eddy viscosity \f$ \nu_\text{t} \f$ is dampened by \f$ f_\mu \f$:
+ * \f[
+ * \nu_\text{t} = C_\mu f_\mu \frac{k^2}{\tilde{\varepsilon}}
+ * \f].
+ *
+ * The auxiliary and dampening functions are defined as:
+ * \f[ D_\varepsilon = 2 \nu \nicefrac{k}{y^2} \f]
+ * \f[ E_\text{k} = -2 \nu \frac{\tilde{\varepsilon}}{y^2} \exp \left( -0.5 y^+ \right) \f]
+ * \f[ f_1 = 1 \f]
+ * \f[ f_2 = 1 - 0.22 \exp \left( - \left( \frac{\mathit{Re}_\text{t}}{6} \right)^2 \right) \f]
+ * \f[ f_\mu = 1 - \exp \left( -0.0115 y^+ \right) \f]
+ * \f[ \mathit{Re}_\text{t} = \frac{k^2}{\nu \tilde{\varepsilon}} \f]
+ * .
+ *
+ * Finally, the model is closed with the following constants:
+ * \f[ \sigma_\text{k} = 1.00 \f]
+ * \f[ \sigma_\varepsilon =1.30 \f]
+ * \f[ C_{1\tilde{\varepsilon}} = 1.35 \f]
+ * \f[ C_{2\tilde{\varepsilon}} = 1.80 \f]
+ * \f[ C_\mu = 0.09 \f]
+ */
+
+#ifndef DUMUX_KEPSILON_MODEL_HH
+#define DUMUX_KEPSILON_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/properties.hh>
+#include <dumux/freeflow/rans/model.hh>
+
+#include "fluxvariables.hh"
+#include "indices.hh"
+#include "localresidual.hh"
+#include "volumevariables.hh"
+#include "vtkoutputfields.hh"
+
+namespace Dumux
+{
+namespace Properties {
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief Traits for the k-epsilon model
+ */
+template<int dimension>
+struct KEpsilonModelTraits : RANSModelTraits<dimension>
+{
+ //! The dimension of the model
+ static constexpr int dim() { return dimension; }
+
+ //! There are as many momentum balance equations as dimensions,
+ //! one mass balance equation and two turbulent transport equations
+ static constexpr int numEq() { return dim()+1+2; }
+
+ //! The number of components
+ static constexpr int numComponents() { return 1; }
+
+ //! the indices
+ using Indices = KEpsilonIndices<dim(), numComponents()>;
+};
+
+///////////////////////////////////////////////////////////////////////////
+// default property values for the isothermal k-epsilon model
+///////////////////////////////////////////////////////////////////////////
+
+//! The type tag for the single-phase, isothermal k-epsilon model
+NEW_TYPE_TAG(KEpsilon, INHERITS_FROM(RANS));
+
+//!< states some specifics of the isothermal k-epsilon model
+SET_PROP(KEpsilon, ModelTraits)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
+ static constexpr int dim = GridView::dimension;
+public:
+ using type = KEpsilonModelTraits<dim>;
+};
+
+//! The flux variables
+SET_PROP(KEpsilon, FluxVariables)
+{
+private:
+ using BaseFluxVariables = NavierStokesFluxVariables<TypeTag>;
+public:
+ using type = KEpsilonFluxVariables<TypeTag, BaseFluxVariables>;
+};
+
+//! The local residual
+SET_PROP(KEpsilon, LocalResidual)
+{
+private:
+ using BaseLocalResidual = NavierStokesResidual<TypeTag>;
+public:
+ using type = KEpsilonResidual<TypeTag, BaseLocalResidual>;
+};
+
+//! Set the volume variables property
+SET_PROP(KEpsilon, VolumeVariables)
+{
+private:
+ using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FSY = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using FST = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+
+ using Traits = NavierStokesVolumeVariablesTraits<PV, FSY, FST, MT>;
+ using NSVolVars = NavierStokesVolumeVariables<Traits>;
+public:
+ using type = KEpsilonVolumeVariables<Traits, NSVolVars>;
+};
+
+//! The specific vtk output fields
+SET_PROP(KEpsilon, VtkOutputFields)
+{
+private:
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+public:
+ using type = KEpsilonVtkOutputFields<FVGridGeometry>;
+};
+
+//////////////////////////////////////////////////////////////////
+// default property values for the non-isothermal k-epsilon model
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the single-phase, isothermal k-epsilon model
+NEW_TYPE_TAG(KEpsilonNI, INHERITS_FROM(RANSNI));
+
+//! The model traits of the non-isothermal model
+SET_PROP(KEpsilonNI, ModelTraits)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
+ static constexpr int dim = GridView::dimension;
+ using IsothermalTraits = KEpsilonModelTraits<dim>;
+public:
+ using type = FreeflowNIModelTraits<IsothermalTraits>;
+};
+
+//! Set the volume variables property
+SET_PROP(KEpsilonNI, VolumeVariables)
+{
+private:
+ using PV = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FSY = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using FST = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+
+ using Traits = NavierStokesVolumeVariablesTraits<PV, FSY, FST, MT>;
+ using NSVolVars = NavierStokesVolumeVariables<Traits>;
+public:
+ using type = KEpsilonVolumeVariables<Traits, NSVolVars>;
+};
+
+//! The specific non-isothermal vtk output fields
+SET_PROP(KEpsilonNI, VtkOutputFields)
+{
+private:
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using IsothermalFields = KEpsilonVtkOutputFields<FVGridGeometry>;
+public:
+ using type = FreeflowNonIsothermalVtkOutputFields<IsothermalFields, ModelTraits>;
+};
+
+// \}
+}
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..db306e02dacf71993ca0d027952978715c763921
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -0,0 +1,134 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup KEpsilonModel
+ * \copydoc Dumux::KEpsilonProblem
+ */
+#ifndef DUMUX_REKEPSILON_PROBLEM_HH
+#define DUMUX_REKEPSILON_PROBLEM_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/staggeredfvproblem.hh>
+#include <dumux/discretization/localview.hh>
+#include <dumux/discretization/staggered/elementsolution.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/rans/problem.hh>
+
+#include "model.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief K-epsilon turbulence problem base class.
+ *
+ * This implements some base functionality for k-epsilon models.
+ */
+template<class TypeTag>
+class KEpsilonProblem : public RANSProblem<TypeTag>
+{
+ using ParentType = RANSProblem<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+
+public:
+ //! The constructor sets the gravity, if desired by the user.
+ KEpsilonProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ useStoredEddyViscosity_ = getParamFromGroup<bool>(this->paramGroup(), "RANS.UseStoredEddyViscosity", true);
+ }
+
+ /*!
+ * \brief Correct size of the static (solution independent) wall variables
+ */
+ void updateStaticWallProperties()
+ {
+ ParentType::updateStaticWallProperties();
+
+ // update size and initial values of the global vectors
+ storedDissipationTilde_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+ storedDynamicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+ storedTurbulentKineticEnergy_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+ }
+
+ /*!
+ * \brief Update the dynamic (solution dependent) relations to the walls
+ *
+ * \param curSol The solution vector.
+ */
+ void updateDynamicWallProperties(const SolutionVector& curSol)
+ {
+ ParentType::updateDynamicWallProperties(curSol);
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const int dofIdx = scv.dofIndex();
+ const auto& cellCenterPriVars = curSol[FVGridGeometry::cellCenterIdx()][dofIdx];
+ PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
+ auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(priVars));
+ // NOTE: first update the turbulence quantities
+ storedDissipationTilde_[elementID] = elemSol[0][Indices::dissipationEqIdx];
+ storedTurbulentKineticEnergy_[elementID] = elemSol[0][Indices::turbulentKineticEnergyEqIdx];
+ // NOTE: then update the volVars
+ VolumeVariables volVars;
+ volVars.update(elemSol, asImp_(), element, scv);
+ storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity();
+ }
+ }
+ }
+
+public:
+ std::vector<Scalar> storedDissipationTilde_;
+ std::vector<Scalar> storedDynamicEddyViscosity_;
+ std::vector<Scalar> storedTurbulentKineticEnergy_;
+ bool useStoredEddyViscosity_;
+
+private:
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation &asImp_()
+ { return *static_cast<Implementation *>(this); }
+
+ //! \copydoc asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation *>(this); }
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/CMakeLists.txt b/dumux/freeflow/rans/twoeq/kepsilon/staggered/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bbbf90a8c6582d728da226f48e959d3d1a164feb
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/CMakeLists.txt
@@ -0,0 +1,5 @@
+#install headers
+install(FILES
+fluxvariables.hh
+localresidual.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/twoeq/kepsilon/staggered)
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..157ab58f3d49992823e77ca6e1cba9943d1df7f0
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
@@ -0,0 +1,172 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+ /*!
+ * \file
+ * \ingroup KEpsilonModel
+ * \copydoc Dumux::KEpsilonFluxVariablesImpl
+ */
+#ifndef DUMUX_KEPSILON_STAGGERED_FLUXVARIABLES_HH
+#define DUMUX_KEPSILON_STAGGERED_FLUXVARIABLES_HH
+
+#include <numeric>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/fluxvariablesbase.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/fluxvariables.hh>
+#include <dumux/freeflow/rans/twoeq/kepsilon/fluxvariables.hh>
+
+namespace Dumux
+{
+
+// forward declaration
+template<class TypeTag, class BaseFluxVariables, DiscretizationMethod discMethod>
+class KEpsilonFluxVariablesImpl;
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief The flux variables class for the k-epsilon model using the staggered grid discretization.
+ */
+template<class TypeTag, class BaseFluxVariables>
+class KEpsilonFluxVariablesImpl<TypeTag, BaseFluxVariables, DiscretizationMethod::staggered>
+: public BaseFluxVariables
+{
+ using ParentType = BaseFluxVariables;
+
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+
+ using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+ using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+ using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+
+ using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
+ using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
+
+ using GridFaceVariables = typename GridVariables::GridFaceVariables;
+ using ElementFaceVariables = typename GridFaceVariables::LocalView;
+ using FaceVariables = typename GridFaceVariables::FaceVariables;
+
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+
+ enum {
+ turbulentKineticEnergyEqIdx = Indices::turbulentKineticEnergyEqIdx,
+ dissipationEqIdx = Indices::dissipationEqIdx,
+ };
+
+public:
+
+ /*!
+ * \brief Computes the flux for the cell center residual.
+ */
+ CellCenterPrimaryVariables computeMassFlux(const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace &scvf,
+ const FluxVariablesCache& fluxVarsCache)
+ {
+ CellCenterPrimaryVariables flux = ParentType::computeMassFlux(problem, element, fvGeometry,
+ elemVolVars, elemFaceVars, scvf, fluxVarsCache);
+
+ // calculate advective flux
+ auto upwindTermK = [](const auto& volVars)
+ {
+ return volVars.turbulentKineticEnergy();
+ };
+ auto upwindTermEpsilon = [](const auto& volVars)
+ {
+ return volVars.dissipationTilde();
+ };
+
+ flux[turbulentKineticEnergyEqIdx - ModelTraits::dim()]
+ = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermK);
+ flux[dissipationEqIdx - ModelTraits::dim()]
+ = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTermEpsilon);
+
+ // calculate diffusive flux
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+ // effective diffusion coefficients
+ Scalar insideCoeff_k = insideVolVars.kinematicViscosity()
+ + insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaK();
+ Scalar outsideCoeff_k = outsideVolVars.kinematicViscosity()
+ + outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaK();
+ Scalar insideCoeff_e = insideVolVars.kinematicViscosity()
+ + insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaEpsilon();
+ Scalar outsideCoeff_e = outsideVolVars.kinematicViscosity()
+ + outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaEpsilon();
+
+ // scale by extrusion factor
+ insideCoeff_k *= insideVolVars.extrusionFactor();
+ outsideCoeff_k *= outsideVolVars.extrusionFactor();
+ insideCoeff_e *= insideVolVars.extrusionFactor();
+ outsideCoeff_e *= outsideVolVars.extrusionFactor();
+
+ // average and distance
+ Scalar coeff_k = (insideCoeff_k + outsideCoeff_k) * 0.5;
+ Scalar coeff_e = (insideCoeff_e + outsideCoeff_e) * 0.5;
+ Scalar distance = 0.0;
+
+ // adapt boundary handling
+ if (scvf.boundary())
+ {
+ distance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
+ coeff_k = insideCoeff_k;
+ coeff_e = insideCoeff_e;
+ }
+ else
+ {
+ distance = (outsideScv.dofPosition() - insideScv.dofPosition()).two_norm();
+ }
+
+ const auto bcTypes = problem.boundaryTypes(element, scvf);
+ if (!(scvf.boundary() && (bcTypes.isOutflow(turbulentKineticEnergyEqIdx)
+ || bcTypes.isSymmetry())))
+ {
+ flux[turbulentKineticEnergyEqIdx - ModelTraits::dim()]
+ += coeff_k / distance
+ * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy())
+ * scvf.area();
+ }
+ if (!(scvf.boundary() && (bcTypes.isOutflow(dissipationEqIdx)
+ || bcTypes.isSymmetry())))
+ {
+ flux[dissipationEqIdx - ModelTraits::dim()]
+ += coeff_e / distance
+ * (insideVolVars.dissipationTilde() - outsideVolVars.dissipationTilde())
+ * scvf.area();
+ }
+ return flux;
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..608a0274b54a28d546ca8f8bb5ca5cf244ae799e
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh
@@ -0,0 +1,128 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+ /*!
+ * \file
+ * \ingroup KepsilonModel
+ * \copydoc Dumux::KepsilonResidualImpl
+ */
+#ifndef DUMUX_STAGGERED_KEPSILON_LOCAL_RESIDUAL_HH
+#define DUMUX_STAGGERED_KEPSILON_LOCAL_RESIDUAL_HH
+
+#include <dune/common/hybridutilities.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/freeflow/navierstokes/localresidual.hh>
+
+namespace Dumux {
+
+// forward declaration
+template<class TypeTag, class BaseLocalResidual, DiscretizationMethod discMethod>
+class KEpsilonResidualImpl;
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief Element-wise calculation of the residual for k-epsilon models using the staggered discretization
+ */
+template<class TypeTag, class BaseLocalResidual>
+class KEpsilonResidualImpl<TypeTag, BaseLocalResidual, DiscretizationMethod::staggered>
+: public BaseLocalResidual
+{
+ using ParentType = BaseLocalResidual;
+
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+
+ using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
+ using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+ using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+
+ using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
+ using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+
+ using GridFaceVariables = typename GridVariables::GridFaceVariables;
+ using ElementFaceVariables = typename GridFaceVariables::LocalView;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+
+public:
+ using ParentType::ParentType;
+
+ //! Evaluate fluxes entering or leaving the cell center control volume.
+ CellCenterPrimaryVariables computeStorageForCellCenter(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
+
+ static constexpr int turbulentKineticEnergyEqIdx = Indices::turbulentKineticEnergyEqIdx - ModelTraits::dim();
+ static constexpr int dissipationEqIdx = Indices::dissipationEqIdx - ModelTraits::dim();
+ storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
+ storage[dissipationEqIdx] = volVars.dissipationTilde();
+
+ return storage;
+ }
+
+ CellCenterPrimaryVariables computeSourceForCellCenter(const Problem& problem,
+ const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolume &scv) const
+ {
+ CellCenterPrimaryVariables source = ParentType::computeSourceForCellCenter(problem, element, fvGeometry,
+ elemVolVars, elemFaceVars, scv);
+
+ const auto& volVars = elemVolVars[scv];
+
+ static constexpr int turbulentKineticEnergyEqIdx = Indices::turbulentKineticEnergyEqIdx - ModelTraits::dim();
+ static constexpr int dissipationEqIdx = Indices::dissipationEqIdx - ModelTraits::dim();
+
+ // production
+ source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.kinematicEddyViscosity()
+ * volVars.stressTensorScalarProduct();
+ source[dissipationEqIdx] += volVars.cOneEpsilon() * volVars.fOne()
+ * volVars.dissipationTilde() / volVars.turbulentKineticEnergy()
+ * 2.0 * volVars.kinematicEddyViscosity()
+ * volVars.stressTensorScalarProduct();
+
+ // destruction
+ source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde();
+ source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo()
+ * volVars.dissipationTilde() * volVars.dissipationTilde()
+ / volVars.turbulentKineticEnergy();
+
+ // dampening functions
+ source[turbulentKineticEnergyEqIdx] -= volVars.dValue();
+ source[dissipationEqIdx] += volVars.eValue();
+
+ return source;
+ }
+};
+}
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/volumevariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..83fa48cc592d7f660a1d9b5175ae63f553663cbc
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/volumevariables.hh
@@ -0,0 +1,293 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup KEpsilonModel
+ *
+ * \copydoc Dumux::KEpsilonVolumeVariables
+ */
+#ifndef DUMUX_KEPSILON_VOLUME_VARIABLES_HH
+#define DUMUX_KEPSILON_VOLUME_VARIABLES_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+#include <dumux/freeflow/rans/volumevariables.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief Volume variables for the isothermal single-phase k-epsilon model.
+ */
+template <class Traits, class NSVolumeVariables>
+class KEpsilonVolumeVariables
+: public RANSVolumeVariables< Traits, KEpsilonVolumeVariables<Traits, NSVolumeVariables> >
+, public NSVolumeVariables
+{
+ using ThisType = KEpsilonVolumeVariables<Traits, NSVolumeVariables>;
+ using RANSParentType = RANSVolumeVariables<Traits, ThisType>;
+ using NavierStokesParentType = NSVolumeVariables;
+
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+
+ static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
+ static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
+
+public:
+ //! export the underlying fluid system
+ using FluidSystem = typename Traits::FluidSystem;
+ //! export the indices type
+ using Indices = typename Traits::ModelTraits::Indices;
+
+ /*!
+ * \brief Update all quantities for a given control volume
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The object specifying the problem which ought to
+ * be simulated
+ * \param element An element which contains part of the control volume
+ * \param scv The sub-control volume
+ */
+ template<class ElementSolution, class Problem, class Element, class SubControlVolume>
+ void update(const ElementSolution &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const SubControlVolume& scv)
+ {
+ NavierStokesParentType::update(elemSol, problem, element, scv);
+ updateRANSProperties(elemSol, problem, element, scv);
+ }
+
+ /*!
+ * \brief Update all turbulent quantities for a given control volume
+ *
+ * Wall and roughness related quantities are stored. Eddy viscosity is set.
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The object specifying the problem which ought to be simulated
+ * \param element An element which contains part of the control volume
+ * \param scv The sub-control volume
+ */
+ template<class ElementSolution, class Problem, class Element, class SubControlVolume>
+ void updateRANSProperties(const ElementSolution &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const SubControlVolume& scv)
+ {
+ RANSParentType::updateRANSProperties(elemSol, problem, element, scv);
+ turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx];
+ dissipationTilde_ = elemSol[0][Indices::dissipationIdx];
+ storedDissipationTilde_ = problem.storedDissipationTilde_[RANSParentType::elementID()];
+ storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()];
+ stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()];
+ if (problem.useStoredEddyViscosity_)
+ dynamicEddyViscosity_ = problem.storedDynamicEddyViscosity_[RANSParentType::elementID()];
+ else
+ dynamicEddyViscosity_ = calculateEddyViscosity();
+ calculateEddyDiffusivity(problem);
+ }
+
+ /*!
+ * \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow
+ */
+ Scalar dynamicEddyViscosity() const
+ { return dynamicEddyViscosity_; }
+
+ /*!
+ * \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
+ * control volume.
+ */
+ Scalar effectiveViscosity() const
+ { return NavierStokesParentType::viscosity() + dynamicEddyViscosity(); }
+
+ /*!
+ * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
+ * of the fluid-flow in the sub-control volume.
+ */
+ template<bool eB = enableEnergyBalance, typename std::enable_if_t<eB, int> = 0>
+ Scalar effectiveThermalConductivity() const
+ {
+ return NavierStokesParentType::thermalConductivity()
+ + RANSParentType::eddyThermalConductivity();
+ }
+
+ /*!
+ * \brief Returns the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$.
+ */
+ Scalar calculateEddyViscosity()
+ {
+ return cMu() * fMu() * turbulentKineticEnergy() * turbulentKineticEnergy()
+ / dissipationTilde() * NavierStokesParentType::density();
+ }
+
+ /*!
+ * \brief Calculates the eddy diffusivity \f$\mathrm{[m^2/s]}\f$ based
+ * on the kinematic eddy viscosity and the turbulent schmidt number
+ */
+ template<class Problem>
+ void calculateEddyDiffusivity(const Problem& problem)
+ {
+ static const auto turbulentSchmidtNumber
+ = getParamFromGroup<Scalar>(problem.paramGroup(),
+ "RANS.TurbulentSchmidtNumber", 1.0);
+ eddyDiffusivity_ = RANSParentType::kinematicEddyViscosity() / turbulentSchmidtNumber;
+ }
+
+ /*!
+ * \brief Returns the turbulent kinetic energy \f$ m^2/s^2 \f$
+ */
+ Scalar turbulentKineticEnergy() const
+ {
+ return turbulentKineticEnergy_;
+ }
+
+ /*!
+ * \brief Returns an effective dissipation \f$ m^2/s^3 \f$
+ */
+ Scalar dissipationTilde() const
+ {
+ return dissipationTilde_;
+ }
+
+ /*!
+ * \brief Returns the turbulent kinetic energy \f$ m^2/s^2 \f$
+ */
+ Scalar storedTurbulentKineticEnergy() const
+ {
+ return storedTurbulentKineticEnergy_;
+ }
+
+ /*!
+ * \brief Returns an effective dissipation \f$ m^2/s^3 \f$
+ */
+ Scalar storedDissipationTilde() const
+ {
+ return storedDissipationTilde_;
+ }
+
+ /*!
+ * \brief Returns the scalar product of the stress tensor
+ */
+ Scalar stressTensorScalarProduct() const
+ {
+ return stressTensorScalarProduct_;
+ }
+
+ //! \brief Returns the \$f Re_\textrm{T} \$f value
+ const Scalar reT() const
+ {
+ return turbulentKineticEnergy() * turbulentKineticEnergy()
+ / RANSParentType::kinematicViscosity() / dissipationTilde();
+ }
+
+ //! \brief Returns the \$f Re_\textrm{y} \$f value
+ const Scalar reY() const
+ {
+ using std::sqrt;
+ return sqrt(turbulentKineticEnergy()) * RANSParentType::wallDistance()
+ / RANSParentType::kinematicViscosity();
+ }
+
+ //! \brief Returns the \$f C_\mu \$f constant
+ const Scalar cMu() const
+ { return 0.09; }
+
+ //! \brief Returns the \$f \sigma_\textrm{k} \$f constant
+ const Scalar sigmaK() const
+ { return 1.0; }
+
+ //! \brief Returns the \$f \sigma_\varepsilon \$f constant
+ const Scalar sigmaEpsilon() const
+ { return 1.3; }
+
+ //! \brief Returns the \$f C_{1\tilde{\varepsilon}} \$f constant
+ const Scalar cOneEpsilon() const
+ { return 1.35; }
+
+ //! \brief Returns the \$f C_{2\tilde{\varepsilon}} \$f constant
+ const Scalar cTwoEpsilon() const
+ { return 1.8; }
+
+ //! \brief Returns the \$f D \$f value
+ const Scalar dValue() const
+ {
+ return 2.0 * RANSParentType::kinematicViscosity() * turbulentKineticEnergy()
+ / RANSParentType::wallDistance() / RANSParentType::wallDistance();
+ }
+
+ //! \brief Returns the \$f f_\mu \$f value
+ const Scalar fMu() const
+ {
+ using std::exp;
+ using std::pow;
+ return 1.0 - exp(-0.0115 * RANSParentType::yPlus());
+ }
+
+ //! \brief Returns the \$f f_1 \$f value
+ const Scalar fOne() const
+ { return 1.0; }
+
+ //! \brief Returns the \$f f_2 \$f value
+ const Scalar fTwo() const
+ {
+ using std::exp;
+ return 1.0 - 0.22 * exp(-1.0 * (reT() * reT() / 6.0 / 6.0));
+ }
+
+ //! \brief Returns the \$f E \$f value
+ const Scalar eValue() const
+ {
+ using std::exp;
+ return -2.0 * RANSParentType::kinematicViscosity() * dissipationTilde()
+ / RANSParentType::wallDistance() / RANSParentType::wallDistance()
+ * exp(-0.5 * RANSParentType::yPlus());
+ }
+
+ /*!
+ * \brief Returns the eddy diffusivity \f$\mathrm{[m^2/s]}\f$
+ */
+ Scalar eddyDiffusivity() const
+ { return eddyDiffusivity_; }
+
+ /*!
+ * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+ *
+ * \param compIIdx the index of the component which diffusive
+ * \param compJIdx the index of the component with respect to which compIIdx diffuses
+ */
+ Scalar effectiveDiffusivity(int compIIdx, int compJIdx = fluidSystemPhaseIdx) const
+ {
+ return NavierStokesParentType::diffusionCoefficient(compIIdx, compJIdx) + eddyDiffusivity();
+ }
+
+protected:
+ Scalar dynamicEddyViscosity_;
+ Scalar eddyDiffusivity_;
+ Scalar turbulentKineticEnergy_;
+ Scalar dissipationTilde_;
+ Scalar storedTurbulentKineticEnergy_;
+ Scalar storedDissipationTilde_;
+ Scalar stressTensorScalarProduct_;
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/vtkoutputfields.hh b/dumux/freeflow/rans/twoeq/kepsilon/vtkoutputfields.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3250528b447e47246fcd2895c19c70d5cb321218
--- /dev/null
+++ b/dumux/freeflow/rans/twoeq/kepsilon/vtkoutputfields.hh
@@ -0,0 +1,61 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup KEpsilonModel
+ * \copydoc Dumux::KEpsilonVtkOutputFields
+ */
+#ifndef DUMUX_KEPSILON_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_KEPSILON_VTK_OUTPUT_FIELDS_HH
+
+#include <dumux/freeflow/rans/vtkoutputfields.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup KEpsilonModel
+ * \brief Adds vtk output fields for the k-epsilon turbulence model
+ */
+template<class FVGridGeometry>
+class KEpsilonVtkOutputFields : public RANSVtkOutputFields<FVGridGeometry>
+{
+ enum { dim = FVGridGeometry::GridView::dimension };
+
+public:
+ //! Initialize the Reynolds-averaged Navier-Stokes specific vtk output fields.
+ template <class VtkOutputModule>
+ static void init(VtkOutputModule& vtk)
+ {
+ RANSVtkOutputFields<FVGridGeometry>::init(vtk);
+ add(vtk);
+ }
+
+ //! Add the KEpsilon specific vtk output fields.
+ template <class VtkOutputModule>
+ static void add(VtkOutputModule& vtk)
+ {
+ vtk.addVolumeVariable([](const auto& v){ return v.turbulentKineticEnergy(); }, "k");
+ vtk.addVolumeVariable([](const auto& v){ return v.dissipationTilde(); }, "epsilon");
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/freeflow/rans/CMakeLists.txt b/test/freeflow/rans/CMakeLists.txt
index 4b69582c8c6def5e2ca12331bdd3a37bc9030087..1be433a716a78c7d57817170c2ff4986a10b1854 100644
--- a/test/freeflow/rans/CMakeLists.txt
+++ b/test/freeflow/rans/CMakeLists.txt
@@ -10,6 +10,16 @@ dune_add_test(NAME test_pipe_laufer_zeroeq
${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_reference-00044.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer_zeroeq test_pipe_laufer_reference.input")
+dune_add_test(NAME test_pipe_laufer_kepsilon
+ SOURCES test_pipe_laufer.cc
+ CMAKE_GUARD HAVE_UMFPACK
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/pipe_laufer_lowrekepsilon.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_reference-00064.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer_kepsilon test_pipe_laufer_reference.input")
+target_compile_definitions(test_pipe_laufer_kepsilon PUBLIC "KEPSILON=1")
+
dune_add_test(NAME test_pipe_laufer_lowrekepsilon
SOURCES test_pipe_laufer.cc
CMAKE_GUARD HAVE_UMFPACK
diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh
index 20d939fe8e8599f572a04ff7c01f055226849369..b718aba8c7847ab97886ba19f7b11558154290f2 100644
--- a/test/freeflow/rans/pipelauferproblem.hh
+++ b/test/freeflow/rans/pipelauferproblem.hh
@@ -27,6 +27,7 @@
#ifndef DUMUX_PIPE_LAUFER_PROBLEM_HH
#define DUMUX_PIPE_LAUFER_PROBLEM_HH
+#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/turbulenceproperties.hh>
#include <dumux/material/fluidsystems/1pgas.hh>
#include <dumux/material/components/air.hh>
@@ -34,11 +35,13 @@
#if LOWREKEPSILON
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh>
#include <dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh>
+#elif KEPSILON
+#include <dumux/freeflow/rans/twoeq/kepsilon/model.hh>
+#include <dumux/freeflow/rans/twoeq/kepsilon/problem.hh>
#else
#include <dumux/freeflow/rans/zeroeq/model.hh>
#include <dumux/freeflow/rans/zeroeq/problem.hh>
#endif
-#include <dumux/discretization/staggered/freeflow/properties.hh>
namespace Dumux
{
@@ -52,6 +55,8 @@ NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNI))
#else
#if LOWREKEPSILON
NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, LowReKEpsilon));
+#elif KEPSILON
+NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, KEpsilon));
#else
NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEq));
#endif
@@ -89,6 +94,10 @@ template <class TypeTag>
class PipeLauferProblem : public LowReKEpsilonProblem<TypeTag>
{
using ParentType = LowReKEpsilonProblem<TypeTag>;
+#elif KEPSILON
+class PipeLauferProblem : public KEpsilonProblem<TypeTag>
+{
+ using ParentType = KEpsilonProblem<TypeTag>;
#else
class PipeLauferProblem : public ZeroEqProblem<TypeTag>
{
@@ -197,7 +206,7 @@ public:
values.setOutflow(Indices::energyBalanceIdx);
#endif
-#if LOWREKEPSILON
+#if LOWREKEPSILON || KEPSILON
values.setOutflow(Indices::turbulentKineticEnergyEqIdx);
values.setOutflow(Indices::dissipationEqIdx);
#endif
@@ -211,7 +220,7 @@ public:
values.setDirichlet(Indices::temperatureIdx);
#endif
-#if LOWREKEPSILON
+#if LOWREKEPSILON || KEPSILON
values.setDirichlet(Indices::turbulentKineticEnergyIdx);
values.setDirichlet(Indices::dissipationIdx);
#endif
@@ -271,7 +280,7 @@ public:
}
#endif
-#if LOWREKEPSILON
+#if LOWREKEPSILON || KEPSILON
if (time() < eps_ && startWithZeroVelocity_)
{
values[Indices::velocityXIdx] = 0.0;