Merge branch 'freeflow/ransncni' into 'master'

Freeflow/ransncni

See merge request !867

dumux/discretization/staggered/freeflow/fourierslaw.hh

@@ -85,8 +85,8 @@ public:
         const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
         // effective conductivity tensors
-        auto insideLambda = insideVolVars.thermalConductivity();
-        auto outsideLambda = outsideVolVars.thermalConductivity();
+        auto insideLambda = insideVolVars.effectiveThermalConductivity();
+        auto outsideLambda = outsideVolVars.effectiveThermalConductivity();
 
         // scale by extrusion factor
         insideLambda *= insideVolVars.extrusionFactor();

dumux/freeflow/navierstokes/volumevariables.hh

@@ -226,7 +226,7 @@ protected:
  */
 template <class TypeTag>
 class NavierStokesVolumeVariablesImplementation<TypeTag, true>
-: public NavierStokesVolumeVariablesImplementation<TypeTag, false>
+: virtual public NavierStokesVolumeVariablesImplementation<TypeTag, false>
 {
     using ParentType = NavierStokesVolumeVariablesImplementation<TypeTag, false>;
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
@@ -326,6 +326,12 @@ public:
     Scalar thermalConductivity() const
     { return FluidSystem::thermalConductivity(this->fluidState_, phaseIdx); }
 
+    /*!
+     * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$.
+     */
+    Scalar effectiveThermalConductivity() const
+    { return thermalConductivity(); }
+
     //! The temperature is a primary variable for non-isothermal models
     using ParentType::temperature;
     template<class ElementSolution>

dumux/freeflow/navierstokesnc/model.hh

@@ -20,20 +20,7 @@
  * \file
  * \ingroup NavierStokesNCModel
  *
- * \brief A single-phase, multi-component isothermal Navier-Stokes model
- *
- * This model implements a single-phase, multi-component isothermal Navier-Stokes model, solving the <B> momentum balance equation </B>
- * \f[
- \frac{\partial (\varrho \textbf{v})}{\partial t} + \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\textup{T}}) = \nabla \cdot (\mu (\nabla \textbf{v} + \nabla \textbf{v}^{\textup{T}}))
-     - \nabla p + \varrho \textbf{g} - \textbf{f}
- * \f]
- * By setting the property <code>EnableInertiaTerms</code> to <code>false</code> the Stokes
- * equation can be solved. In this case the term
- * \f[
- *    \nabla \cdot (\varrho \textbf{v} \textbf{v}^{\textup{T}})
- * \f]
- * is neglected.
- *
+ * \copydoc Dumux::NavierStokesModel
  *
  * The system is closed by a <B> component mass/mole balance equation </B> for each component \f$\kappa\f$:
  * \f[

dumux/freeflow/nonisothermal/model.hh

@@ -26,9 +26,18 @@
  * \f[
  *    \frac{\partial (\varrho  v)}{\partial t}
  *    + \nabla \cdot \left( \varrho h {\boldsymbol{v}}
- *    - \lambda \textbf{grad}\, T \right) - q_T = 0
+ *    - \lambda_\text{eff} \textbf{grad}\, T \right) - q_T = 0
  * \f]
  *
+ *
+ * For laminar Navier-Stokes flow the effective thermal conductivity is the fluid
+ * thermal conductivity: \f$ \lambda_\text{eff} = \lambda \f$.
+ *
+ * For turbulent Reynolds-averaged Navier-Stokes flow the eddy thermal conductivity is added:
+ *  \f$ \lambda_\text{eff} = \lambda + \lambda_\text{t} \f$.
+ * The eddy thermal conductivity \f$ \lambda_\text{t} \f$ is related to the eddy viscosity \f$ \nu_\text{t} \f$
+ * by the turbulent Prandtl number:
+ * \f[ \lambda_\text{t} = \frac{\nu_\text{t} \varrho c_\text{p}}{\mathrm{Pr}_\text{t}} \f]
  */
 
 #ifndef DUMUX_STAGGERED_NI_MODEL_HH

dumux/freeflow/nonisothermal/ransvtkoutputfields.hh

+// -*- 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 RANSNIModel
+ * \copydoc Dumux::RANSNonIsothermalVtkOutputFields
+ */
+#ifndef DUMUX_RANS_NI_OUTPUT_FIELDS_HH
+#define DUMUX_RANS_NI_OUTPUT_FIELDS_HH
+
+#include <dumux/common/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup RANSNIModel
+ * \brief Adds vtk output fields specific to non-isothermal RANS models
+ */
+template<class NavierStokesVtkOutputFields, class RANSVtkOutputFields>
+class RANSNonIsothermalVtkOutputFields
+{
+
+public:
+    //! Initialize the non-isothermal specific vtk output fields.
+    template <class VtkOutputModule>
+    static void init(VtkOutputModule& vtk)
+    {
+        NavierStokesVtkOutputFields::init(vtk);
+        add(vtk);
+    }
+
+    //! Add the RANS specific vtk output fields.
+    template <class VtkOutputModule>
+    static void add(VtkOutputModule& vtk)
+    {
+        RANSVtkOutputFields::add(vtk);
+        vtk.addVolumeVariable( [](const auto& v){ return v.temperature(); }, "temperature");
+        vtk.addVolumeVariable( [](const auto& v){ return v.eddyThermalConductivity(); }, "lambda_t");
+    }
+};
+
+} // end namespace Dumux
+
+#endif

dumux/freeflow/rans/model.hh

@@ -41,7 +41,9 @@
 #include <dumux/discretization/methods.hh>
 #include <dumux/freeflow/properties.hh>
 #include <dumux/freeflow/navierstokes/model.hh>
-#include <dumux/freeflow/nonisothermal/model.hh>
+#include <dumux/freeflow/navierstokes/indices.hh>
+#include <dumux/freeflow/nonisothermal/indices.hh>
+#include <dumux/freeflow/nonisothermal/ransvtkoutputfields.hh>
 #include <dumux/material/fluidstates/immiscible.hh>
 
 #include "volumevariables.hh"
@@ -77,8 +79,52 @@ SET_PROP(RANS, VtkOutputFields)
 private:
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
 public:
-     using type = RANSVtkOutputFields<FVGridGeometry>;
+    using type = RANSVtkOutputFields<FVGridGeometry>;
 };
+
+//////////////////////////////////////////////////////////////////
+// Property values for non-isothermal Reynolds-averaged Navier-Stokes model
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the single-phase, isothermal Reynolds-Averaged Navier-Stokes model
+NEW_TYPE_TAG(RANSNI, INHERITS_FROM(RANS));
+
+//! The model traits of the non-isothermal model
+SET_PROP(RANSNI, ModelTraits)
+{
+private:
+    using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
+    static constexpr int dim = GridView::dimension;
+    using IsothermalTraits = NavierStokesModelTraits<dim>;
+public:
+    using type = NavierStokesNIModelTraits<IsothermalTraits>;
+};
+
+//! The indices required by the non-isothermal single-phase model
+SET_PROP(RANSNI, Indices)
+{
+private:
+    static constexpr int numEq = GET_PROP_TYPE(TypeTag, ModelTraits)::numEq();
+    static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
+    using IsothermalIndices = NavierStokesIndices<dim, numEq>;
+public:
+    using type = NavierStokesNonIsothermalIndices<dim, numEq, IsothermalIndices>;
+};
+
+//! The specific non-isothermal vtk output fields
+SET_PROP(RANSNI, VtkOutputFields)
+{
+private:
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using NavierStokesFields = NavierStokesVtkOutputFields<FVGridGeometry>;
+    using RANSFields = RANSVtkOutputFields<FVGridGeometry>;
+public:
+    using type = RANSNonIsothermalVtkOutputFields<NavierStokesFields, RANSFields>;
+};
+
+//! Use Fourier's Law as default heat conduction type
+SET_TYPE_PROP(RANSNI, HeatConductionType, FouriersLaw<TypeTag>);
+
 // \}
 }
 

dumux/freeflow/rans/volumevariables.hh

@@ -29,6 +29,7 @@
 #include <dune/common/fmatrix.hh>
 
 #include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
 #include <dumux/material/fluidstates/immiscible.hh>
 #include <dumux/freeflow/navierstokes/volumevariables.hh>
 
@@ -71,31 +72,26 @@ class RANSVolumeVariablesImplementation<TypeTag, false>
     using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 
 public:
-
-    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
-
     /*!
-     * \brief Update all quantities for a given control volume
+     * \brief Update all turbulent quantities for a given control volume
      *
      * Wall related quantities are stored and the calculateEddyViscosity(...)
      * function of the turbulence model implementation is called.
      *
      * \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 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>
-    void update(const ElementSolution &elemSol,
-                const Problem &problem,
-                const Element &element,
-                const SubControlVolume& scv)
+    void updateRANSProperties(const ElementSolution &elemSol,
+                              const Problem &problem,
+                              const Element &element,
+                              const SubControlVolume& scv)
     {
         using std::abs;
         using std::max;
         using std::sqrt;
-        ParentType::update(elemSol, problem, element, scv);
 
         // calculate characteristic properties of the turbulent flow
         elementID_ = problem.fvGridGeometry().elementMapper().index(element);
@@ -109,9 +105,10 @@ public:
         uStar_ = sqrt(problem.kinematicViscosity_[wallElementID_]
                       * abs(problem.velocityGradients_[wallElementID_][flowNormalAxis][wallNormalAxis]));
         uStar_ = max(uStar_, 1e-10); // zero values lead to numerical problems in some turbulence models
-        yPlus_ = wallDistance_ * uStar_ / asImp_().kinematicViscosity();
+        yPlus_ = wallDistance_ * uStar_ / kinematicViscosity();
         uPlus_ = velocity_[flowNormalAxis] / uStar_;
-    };
+        dynamicEddyViscosity_ = 0.0; // will be set by the specific RANS implementation
+    }
 
     /*!
      * \brief Return the element ID of the control volume.
@@ -166,9 +163,15 @@ public:
      *        control volume.
      */
     Scalar dynamicEddyViscosity() const
+    { return dynamicEddyViscosity_; }
+
+    /*!
+     * \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
+     *        control volume.
+     */
+    void setDynamicEddyViscosity(Scalar value)
     {
-        DUNE_THROW(Dune::NotImplemented,
-                   "The dynamicEddyViscosity() has to specified by the RANS implementation.");
+        dynamicEddyViscosity_ = value;
     }
 
     /*!
@@ -176,30 +179,21 @@ public:
      *        control volume.
      */
     Scalar effectiveViscosity() const
-    { return asImp_().viscosity() + asImp_().dynamicEddyViscosity(); }
+    { return ParentType::viscosity() + dynamicEddyViscosity(); }
 
     /*!
      * \brief Return the kinematic viscosity \f$\mathrm{[m^2/s]}\f$ of the fluid within the
      *        control volume.
      */
     Scalar kinematicViscosity() const
-    { return asImp_().viscosity() / asImp_().density(); }
+    { return ParentType::viscosity() / ParentType::density(); }
 
     /*!
      * \brief Return the kinematic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
      *        control volume.
      */
     Scalar kinematicEddyViscosity() const
-    { return asImp_().dynamicEddyViscosity() / asImp_().density(); }
-
-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); }
+    { return dynamicEddyViscosity() / ParentType::density(); }
 
 protected:
     DimVector velocity_;
@@ -220,49 +214,38 @@ protected:
  */
 template <class TypeTag>
 class RANSVolumeVariablesImplementation<TypeTag, true>
-: public NavierStokesVolumeVariablesImplementation<TypeTag, true>,
-  public RANSVolumeVariablesImplementation<TypeTag, false>
+: virtual public NavierStokesVolumeVariablesImplementation<TypeTag, true>,
+  virtual public RANSVolumeVariablesImplementation<TypeTag, false>
 {
     using ParentTypeNonIsothermal = NavierStokesVolumeVariablesImplementation<TypeTag, true>;
     using ParentTypeIsothermal = RANSVolumeVariablesImplementation<TypeTag, false>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Element = typename GridView::template Codim<0>::Entity;
-    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-
-//     static const int temperatureIdx = Indices::temperatureIdx;
 
 public:
-
-    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
-
     /*!
-     * \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
+     * \brief Calculates the eddy thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ based
+     *        on the kinematic eddy viscosity and the turbulent prandtl number
      */
-    template<class ElementSolution>
-    void update(const ElementSolution &elemSol,
-                const Problem &problem,
-                const Element &element,
-                const SubControlVolume &scv)
+    void calculateEddyThermalConductivity()
     {
-        ParentTypeIsothermal::update(elemSol, problem, element, scv);
-        ParentTypeNonIsothermal::update(elemSol, problem, element, scv);
-        // TODO: convert eddyViscosity to eddyThermalConductivity
+        static const auto turbulentPrandtlNumber
+            = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup),
+                                        "RANS.TurbulentPrandtlNumber", 1.0);
+        eddyThermalConductivity_ = ParentTypeIsothermal::kinematicEddyViscosity()
+                                   * ParentTypeIsothermal::density()
+                                   * ParentTypeNonIsothermal::heatCapacity()
+                                   / turbulentPrandtlNumber;
     }
 
     /*!
      * \brief Returns the eddy thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
-     *        of the flow phase in the sub-control volume.
+     *        of the flow in the sub-control volume.
      */
     Scalar eddyThermalConductivity() const
     { return eddyThermalConductivity_; }
@@ -273,7 +256,7 @@ public:
      */
     Scalar effectiveThermalConductivity() const
     {
-        return FluidSystem::thermalConductivity(this->fluidState_)
+        return ParentTypeNonIsothermal::thermalConductivity()
                + eddyThermalConductivity();
     }
 

dumux/freeflow/rans/vtkoutputfields.hh

@@ -47,7 +47,7 @@ public:
         add(vtk);
     }
 
-    //! Initialize the Navier-Stokes specific vtk output fields.
+    //! Add the RANS specific vtk output fields.
     template <class VtkOutputModule>
     static void add(VtkOutputModule& vtk)
     {

dumux/freeflow/rans/zeroeq/model.hh

@@ -41,27 +41,28 @@
 
 namespace Dumux
 {
+namespace Properties {
 
-// \{
 ///////////////////////////////////////////////////////////////////////////
-// properties for the single-phase Reynolds-Averaged Navier-Stokes 0-Eq. model
+// default property values for the isothermal RANS 0-Eq. model
 ///////////////////////////////////////////////////////////////////////////
-namespace Properties {
-
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
 
 //! The type tag for the single-phase, isothermal Reynolds-Averaged Navier-Stokes 0-Eq. model
 NEW_TYPE_TAG(ZeroEq, INHERITS_FROM(RANS));
 
-///////////////////////////////////////////////////////////////////////////
-// default property values for the isothermal single phase model
-///////////////////////////////////////////////////////////////////////////
-
 //! The volume variables
 SET_TYPE_PROP(ZeroEq, VolumeVariables, ZeroEqVolumeVariables<TypeTag>);
 
+//////////////////////////////////////////////////////////////////
+// default property values for the non-isothermal RANS 0-Eq. model
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the single-phase, isothermal Reynolds-Averaged Navier-Stokes model
+NEW_TYPE_TAG(ZeroEqNI, INHERITS_FROM(RANSNI));
+
+//! The volume variables
+SET_TYPE_PROP(ZeroEqNI, VolumeVariables, ZeroEqVolumeVariables<TypeTag>);
+
 // \}
 }
 

dumux/freeflow/rans/zeroeq/volumevariables.hh

@@ -53,7 +53,7 @@ using ZeroEqVolumeVariables = ZeroEqVolumeVariablesImplementation<TypeTag, GET_P
  */
 template <class TypeTag>
 class ZeroEqVolumeVariablesImplementation<TypeTag, false>
-: public RANSVolumeVariablesImplementation<TypeTag, false>
+: virtual public RANSVolumeVariablesImplementation<TypeTag, false>
 {
     using ParentType = RANSVolumeVariablesImplementation<TypeTag, false>;
     using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
@@ -85,25 +85,44 @@ public:
                 const SubControlVolume& scv)
     {
         ParentType::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>
+    void updateRANSProperties(const ElementSolution &elemSol,
+                              const Problem &problem,
+                              const Element &element,
+                              const SubControlVolume& scv)
+    {
+        ParentType::updateRANSProperties(elemSol, problem, element, scv);
         additionalRoughnessLength_ = problem.additionalRoughnessLength_[this->elementID()];
         yPlusRough_ = wallDistanceRough() * this->uStar() / this->kinematicViscosity();
-        dynamicEddyViscosity_ = calculateEddyViscosity(elemSol, problem, element, scv);
+        calculateEddyViscosity(elemSol, problem, element, scv);
     }
 
     /*!
      * \brief Calculate and set the dynamic eddy viscosity.
      *
      * \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 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>
-    Scalar calculateEddyViscosity(const ElementSolution &elemSol,
-                                  const Problem &problem,
-                                  const Element &element,
-                                  const SubControlVolume& scv)
+    void calculateEddyViscosity(const ElementSolution &elemSol,
+                                const Problem &problem,
+                                const Element &element,
+                                const SubControlVolume& scv)
     {
         using std::abs;
         using std::exp;
@@ -114,7 +133,7 @@ public:
         unsigned int elementID = problem.fvGridGeometry().elementMapper().index(element);
         unsigned int flowNormalAxis = problem.flowNormalAxis_[elementID];
         unsigned int wallNormalAxis = problem.wallNormalAxis_[elementID];
-        Scalar velGrad = abs(asImp_().velocityGradients()[flowNormalAxis][wallNormalAxis]);
+        Scalar velGrad = abs(ParentType::velocityGradients()[flowNormalAxis][wallNormalAxis]);
 
         if (eddyViscosityModel == EddyViscosityModels::none)
         {
@@ -122,14 +141,14 @@ public:
         }
         else if (eddyViscosityModel == EddyViscosityModels::prandtl)
         {
-            Scalar mixingLength = problem.karmanConstant() * asImp_().wallDistanceRough();
+            Scalar mixingLength = problem.karmanConstant() * wallDistanceRough();
             kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
         }
         else if (eddyViscosityModel == EddyViscosityModels::modifiedVanDriest)
         {
-            Scalar mixingLength = problem.karmanConstant() * asImp_().wallDistanceRough()
-                                  * (1.0 - exp(-asImp_().yPlusRough() / 26.0))
-                                  / sqrt(1.0 - exp(-0.26 * asImp_().yPlusRough()));
+            Scalar mixingLength = problem.karmanConstant() * wallDistanceRough()
+                                  * (1.0 - exp(-yPlusRough() / 26.0))
+                                  / sqrt(1.0 - exp(-0.26 * yPlusRough()));
             kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
         }
         else if (eddyViscosityModel == EddyViscosityModels::baldwinLomax)
@@ -141,7 +160,7 @@ public:
             DUNE_THROW(Dune::NotImplemented,
                        "This eddy viscosity model is not implemented: " << eddyViscosityModel);
         }
-        return kinematicEddyViscosity * asImp_().density();
+        ParentType::setDynamicEddyViscosity(kinematicEddyViscosity * ParentType::density());
     }
 
     /*!
@@ -156,26 +175,9 @@ public:
     Scalar yPlusRough() const
     { return yPlusRough_; }
 
-    /*!
-     * \brief Return the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ of the flow within the
-     *        control volume.
-     */
-    Scalar dynamicEddyViscosity() const
-    { return dynamicEddyViscosity_; }
-
-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); }
-
 protected:
     Scalar additionalRoughnessLength_;
     Scalar yPlusRough_;
-    Scalar dynamicEddyViscosity_;
 };
 
 /*!
@@ -184,8 +186,58 @@ protected:
  */
 template <class TypeTag>
 class ZeroEqVolumeVariablesImplementation<TypeTag, true>
-: public RANSVolumeVariablesImplementation<TypeTag, true>
-{ };
+: virtual public ZeroEqVolumeVariablesImplementation<TypeTag, false>,
+  virtual public RANSVolumeVariablesImplementation<TypeTag, true>
+{
+    using ParentTypeNonIsothermal = RANSVolumeVariablesImplementation<TypeTag, true>;
+    using ParentTypeIsothermal = ZeroEqVolumeVariablesImplementation<TypeTag, false>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+    /*!
+     * \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>
+    void update(const ElementSolution &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume &scv)
+    {
+        ParentTypeNonIsothermal::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>
+    void updateRANSProperties(const ElementSolution &elemSol,
+                              const Problem &problem,
+                              const Element &element,
+                              const SubControlVolume& scv)
+    {
+        ParentTypeIsothermal::updateRANSProperties(elemSol, problem, element, scv);
+        ParentTypeNonIsothermal::calculateEddyThermalConductivity();
+    }
+};
 }
 
 #endif

dumux/freeflow/ransnc/model.hh

@@ -22,7 +22,7 @@
  *
  * \brief A single-phase, multi-component isothermal Navier-Stokes model
  *
- * \copydoc RANSModel
+ * \copydoc Dumux::RANSModel
  *
  * The system is closed by a <B> component mass/mole balance equation </B> for each component \f$\kappa\f$:
  * \f[
@@ -69,14 +69,14 @@ namespace Dumux {
 
 /*!
  * \ingroup RANSModel
- * \brief Traits for the Reynolds-averaged Navier-Stokes multi-component model
+ * \brief Traits for the RANS multi-component model
  */
 template<int dim, int nComp>
 struct RANSNCModelTraits : NavierStokesNCModelTraits<dim, nComp>
 { };
 
 ///////////////////////////////////////////////////////////////////////////
-// properties for the single-phase, multi-component Reynolds-averaged Navier-Stokes model
+// properties for the single-phase, multi-component RANS model
 ///////////////////////////////////////////////////////////////////////////
 namespace Properties {
 
@@ -84,13 +84,10 @@ namespace Properties {
 // Type tags
 //////////////////////////////////////////////////////////////////
 
-//! The type tag for the single-phase, multi-component isothermal Reynolds-averaged Navier-Stokes model
+//! The type tags for the single-phase, multi-component isothermal RANS model
 NEW_TYPE_TAG(RANSNC, INHERITS_FROM(RANS, NavierStokesNC));
 NEW_TYPE_TAG(ZeroEqNC, INHERITS_FROM(ZeroEq, RANSNC));
 
-// //! The type tag for the single-phase, multi-component non-isothermal Reynolds-averaged Navier-Stokes model
-// NEW_TYPE_TAG(RANSNCNI, INHERITS_FROM(RANS, NavierStokesNC));
-
 ///////////////////////////////////////////////////////////////////////////
 // default property values
 ///////////////////////////////////////////////////////////////////////////
@@ -124,23 +121,56 @@ public:
 };
 
 //////////////////////////////////////////////////////////////////////////
-// Property values for non-isothermal multi-component Reynolds-averaged Navier-Stokes model
+// Property values for non-isothermal multi-component RANS model
 //////////////////////////////////////////////////////////////////////////
 
-// //! The non-isothermal vtk output fields
-// SET_PROP(RANSNCNI, VtkOutputFields)
-// {
-// private:
-//     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
-//     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-//     static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-//     using IsothermalFields = RANSNCVtkOutputFields<FVGridGeometry, FluidSystem, phaseIdx>;
-// public:
-//      using type = NavierStokesNonIsothermalVtkOutputFields<IsothermalFields>;
-// };
-//
-// //! Use Fourier's Law as default heat conduction type
-// SET_TYPE_PROP(RANSNCNI, HeatConductionType, FouriersLaw<TypeTag>);
+//! The type tags for the single-phase, multi-component non-isothermal RANS models
+NEW_TYPE_TAG(RANSNCNI, INHERITS_FROM(RANSNI, NavierStokesNC));
+NEW_TYPE_TAG(ZeroEqNCNI, INHERITS_FROM(ZeroEqNI, RANSNCNI));
+
+//! The volume variables
+SET_TYPE_PROP(RANSNCNI, VolumeVariables, RANSNCVolumeVariables<TypeTag>);
+SET_TYPE_PROP(ZeroEqNCNI, SinglePhaseVolumeVariables, ZeroEqVolumeVariables<TypeTag>);
+
+//! The model traits of the non-isothermal model
+SET_PROP(RANSNCNI, ModelTraits)
+{
+private:
+    using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
+    static constexpr int dim = GridView::dimension;
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static constexpr int numComponents = FluidSystem::numComponents;
+    using IsothermalModelTraits = NavierStokesNCModelTraits<dim, numComponents>;
+public:
+    using type = NavierStokesNIModelTraits<IsothermalModelTraits>;
+};
+
+//! The indices
+SET_PROP(RANSNCNI, Indices)
+{
+private:
+    static constexpr int numEq = GET_PROP_TYPE(TypeTag, ModelTraits)::numEq();
+    static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    using IsothermalIndices = NavierStokesNCIndices<dim, numEq, phaseIdx, replaceCompEqIdx>;
+public:
+    using type = NavierStokesNonIsothermalIndices<dim, numEq, IsothermalIndices>;
+};
+
+//! The specific vtk output fields
+SET_PROP(ZeroEqNCNI, VtkOutputFields)
+{
+private:
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using NavierStokesFields = NavierStokesVtkOutputFields<FVGridGeometry>;
+    using RANSFields = RANSVtkOutputFields<FVGridGeometry>;
+    using SinglePhaseVtkOutputFields = RANSNonIsothermalVtkOutputFields<NavierStokesFields, RANSFields>;
+public:
+    using type = RANSNCVtkOutputFields<FVGridGeometry, FluidSystem, phaseIdx, SinglePhaseVtkOutputFields>;
+};
 
 // \}
 } // end namespace Properties

dumux/freeflow/ransnc/volumevariables.hh

@@ -38,8 +38,8 @@ namespace Dumux {
  */
 template <class TypeTag>
 class RANSNCVolumeVariables
-      : public GET_PROP_TYPE(TypeTag, SinglePhaseVolumeVariables),
-        public NavierStokesNCVolumeVariables<TypeTag>
+    : virtual public GET_PROP_TYPE(TypeTag, SinglePhaseVolumeVariables),
+      virtual public NavierStokesNCVolumeVariables<TypeTag>
 {
     using ParentTypeSinglePhase = typename GET_PROP_TYPE(TypeTag, SinglePhaseVolumeVariables);
     using ParentTypeCompositional = NavierStokesNCVolumeVariables<TypeTag>;
@@ -54,7 +54,6 @@ class RANSNCVolumeVariables
     static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 
 public:
-
     /*!
      * \brief Update all quantities for a given control volume
      *
@@ -70,12 +69,17 @@ public:
                 const Element &element,
                 const SubControlVolume& scv)
     {
-        this->extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol);
-
-        // NOTE: the ParentTypeCompositional has to be called first
         ParentTypeCompositional::update(elemSol, problem, element, scv);
-        ParentTypeSinglePhase::update(elemSol, problem, element, scv);
+        ParentTypeSinglePhase::updateRANSProperties(elemSol, problem, element, scv);
+        calculateEddyDiffusivity();
+    }
 
+    /*!
+     * \brief Calculates the eddy diffusivity \f$\mathrm{[m^2/s]}\f$ based
+     *        on the kinematic eddy viscosity and the turbulent schmidt number
+     */
+    void calculateEddyDiffusivity()
+    {
         static const auto turbulentSchmidtNumber
             = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup),
                                         "RANS.TurbulentSchmidtNumber", 1.0);
@@ -103,13 +107,6 @@ public:
     }
 
 protected:
-
-    Implementation &asImp_()
-    { return *static_cast<Implementation*>(this); }
-
-    const Implementation &asImp_() const
-    { return *static_cast<const Implementation*>(this); }
-
     Scalar eddyDiffusivity_;
 };
 

test/freeflow/rans/CMakeLists.txt

 add_input_file_links()
 dune_symlink_to_source_files(FILES laufer_re50000.csv laufer_re50000_u+y+.csv)
 
-dune_add_test(NAME test_pipe_laufer
+dune_add_test(NAME test_pipe_laufer_zeroeq
               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_zeroeq.vtu
-                                ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_zeroeq_reference-00021.vtu
-                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer test_pipe_laufer_reference.input")
+                                ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer_zeroeq_reference-00020.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer_zeroeq test_pipe_laufer_reference.input")
+
+dune_add_test(NAME test_pipe_zeroeqni
+              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_zeroeqni.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/pipe_zeroeqni-00034.vtu
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_zeroeqni test_pipe_zeroeqni.input")
+target_compile_definitions(test_pipe_zeroeqni PUBLIC "NONISOTHERMAL=1")

test/freeflow/rans/pipelauferproblem.hh

@@ -41,7 +41,11 @@ class PipeLauferProblem;
 
 namespace Properties
 {
+#if NONISOTHERMAL
+NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNI));
+#else
 NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEq));
+#endif
 
 // the fluid system
 SET_PROP(PipeLauferProblem, FluidSystem)
@@ -82,6 +86,10 @@ class PipeLauferProblem : public ZeroEqProblem<TypeTag>
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     enum { dimWorld = GridView::dimensionworld };
     enum {
+#if NONISOTHERMAL
+        temperatureIdx = Indices::temperatureIdx,
+        energyBalanceIdx = Indices::energyBalanceIdx,
+#endif
         totalMassBalanceIdx = Indices::totalMassBalanceIdx,
         momentumBalanceIdx = Indices::momentumBalanceIdx,
         pressureIdx = Indices::pressureIdx,
@@ -105,6 +113,8 @@ public:
     : ParentType(fvGridGeometry), eps_(1e-6)
     {
         inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
+        inletTemperature_ = getParam<Scalar>("Problem.InletTemperature", 283.15);
+        wallTemperature_ = getParam<Scalar>("Problem.WallTemperature", 323.15);
         sandGrainRoughness_ = getParam<Scalar>("Problem.SandGrainRoughness", 0.0);
     }
 
@@ -115,9 +125,8 @@ public:
 
     bool isOnWall(const GlobalPosition &globalPos) const
     {
-        Scalar localEps_ = 1e-6; // cannot use the epsilon, because ParentType is initialized first
-        return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + localEps_
-              || globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - localEps_;
+        return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_
+               || globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
     }
 
     Scalar sandGrainRoughnessAtPos(const GlobalPosition &globalPos) const
@@ -131,12 +140,10 @@ public:
     }
 
    /*!
-     * \brief Return the temperature within the domain in [K].
-     *
-     * This problem assumes a temperature of 10 degrees Celsius.
+     * \brief Return the temperature [K] within the domain for the isothermal model.
      */
     Scalar temperature() const
-    { return 273.15 + 10; } // 10C
+    { return inletTemperature_; }
 
    /*!
      * \brief Return the sources within the domain.
@@ -165,6 +172,13 @@ public:
 
         // set Dirichlet values for the velocity everywhere
         values.setDirichlet(momentumBalanceIdx);
+#if NONISOTHERMAL
+        values.setDirichlet(energyBalanceIdx);
+        if (isOutlet(globalPos))
+        {
+            values.setOutflow(energyBalanceIdx);
+        }
+#endif
 
         // set a fixed pressure in one cell
         if (isOutlet(globalPos))
@@ -186,7 +200,18 @@ public:
      */
     PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
     {
-        return initialAtPos(globalPos);
+        PrimaryVariables values(initialAtPos(globalPos));
+#if NONISOTHERMAL
+        if (time() > 10.0)
+        {
+            values[temperatureIdx] = inletTemperature_;
+            if (isOnWall(globalPos))
+            {
+                values[temperatureIdx] = wallTemperature_;
+            }
+        }
+#endif
+        return values;
     }
 
     // \}
@@ -205,9 +230,17 @@ public:
     {
         PrimaryVariables values(0.0);
         values[pressureIdx] = 1.0e+5;
-        if(!isOnWall(globalPos))
+        values[velocityXIdx] = inletVelocity_;
+#if NONISOTHERMAL
+        values[temperatureIdx] = inletTemperature_;
+        if (isOnWall(globalPos))
+        {
+            values[temperatureIdx] = wallTemperature_;
+        }
+#endif
+        if (isOnWall(globalPos))
         {
-            values[velocityXIdx] = inletVelocity_;
+            values[velocityXIdx] = 0.0;
         }
 
         return values;
@@ -217,8 +250,6 @@ public:
     void setTimeLoop(TimeLoopPtr timeLoop)
     {
         timeLoop_ = timeLoop;
-        if(inletVelocity_ > eps_)
-            timeLoop_->setCheckPoint({10.0, 50.0, 100.});
     }
 
     Scalar time() const
@@ -227,6 +258,11 @@ public:
     }
 
 private:
+    bool isInlet(const GlobalPosition& globalPos) const
+    {
+        return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_;
+    }
+
     bool isOutlet(const GlobalPosition& globalPos) const
     {
         return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
@@ -234,6 +270,8 @@ private:
 
     Scalar eps_;
     Scalar inletVelocity_;
+    Scalar inletTemperature_;
+    Scalar wallTemperature_;
     Scalar sandGrainRoughness_;
     TimeLoopPtr timeLoop_;
 };

test/freeflow/rans/test_pipe_zeroeqni.input

+[TimeLoop]
+DtInitial = 10e-3 # [s]
+TEnd = 43200 # [s]
+
+[Grid]
+Verbosity = true
+Positions0 = 0.0 10.0
+Positions1 = 0.0 1.2345 2.469
+Cells0 = 16
+Cells1 = 8 8
+Grading1 = 1.5 -1.5
+
+[Problem]
+Name = pipe_zeroeqni
+InletVelocity = 0.25 # [m/s]
+InletTemperature = 283.15 # [K]
+WallTemperature = 303.15 # [K]
+EnableGravity = false
+
+[RANS]
+EddyViscosityModel = 3
+
+[Newton]
+MaxSteps = 10
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false

test/freeflow/ransnc/CMakeLists.txt

@@ -6,5 +6,15 @@ dune_add_test(NAME test_channel_zeroeq2c
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_channel_zeroeq2c.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2c-00033.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2c-00031.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2c test_channel_zeroeq2c.input")
+
+dune_add_test(NAME test_channel_zeroeq2cni
+              SOURCES test_channel_zeroeq.cc
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS       --script fuzzy
+                             --files ${CMAKE_SOURCE_DIR}/test/references/test_channel_zeroeq2cni.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2cni-00018.vtu
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2cni test_channel_zeroeq2cni.input")
+target_compile_definitions(test_channel_zeroeq2cni PUBLIC "NONISOTHERMAL=1")

test/freeflow/ransnc/channeltestproblem.hh

@@ -39,10 +39,10 @@ class ChannelNCTestProblem;
 namespace Properties
 {
 
-#if !NONISOTHERMAL
-NEW_TYPE_TAG(ChannelNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNC));
-#else
+#if NONISOTHERMAL
 NEW_TYPE_TAG(ChannelNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNCNI));
+#else
+NEW_TYPE_TAG(ChannelNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNC));
 #endif
 
 NEW_PROP_TAG(FluidSystem);
@@ -203,7 +203,7 @@ public:
             values.setOutflow(totalMassBalanceIdx);
             values.setOutflow(transportEqIdx);
 #if NONISOTHERMAL
-            values.setOutflow(energyBalanceIdx);
+            values.setDirichlet(energyBalanceIdx);
 #endif
         }
         else if (isSymmetry(globalPos))
@@ -223,18 +223,21 @@ public:
     {
         PrimaryVariables values = initialAtPos(globalPos);
 
-        if (isInlet(globalPos)
-            && globalPos[1] > 0.4 * this->fvGridGeometry().bBoxMax()[1]
-            && globalPos[1] < 0.6 * this->fvGridGeometry().bBoxMax()[1])
+        if (time() > 10.0)
         {
-            values[transportCompIdx] = 1e-3;
-        }
+            if (isInlet(globalPos)
+                && globalPos[1] > 0.4 * this->fvGridGeometry().bBoxMax()[1]
+                && globalPos[1] < 0.6 * this->fvGridGeometry().bBoxMax()[1])
+            {
+                values[transportCompIdx] = 1e-3;
+            }
 #if NONISOTHERMAL
-        if (isOnWall(globalPos))
-        {
-            values[temperatureIdx] += 30.;
-        }
+            if (isOnWall(globalPos))
+            {
+                values[temperatureIdx] += 30.;
+            }
 #endif
+        }
 
         return values;
     }

test/freeflow/ransnc/test_channel_zeroeq2cni.input

+[TimeLoop]
+DtInitial = 1 # [s]
+TEnd = 1000 # [s]
+
+[Grid]
+Verbosity = true
+Positions0 = 0.0 6.0
+Positions1 = 0.0 0.5
+Cells0 = 10
+Cells1 = 20
+Grading1 = 1.4
+
+[Problem]
+Name = test_channel_zeroeq2cni # name passed to the output routines
+InletVelocity = 0.1 # [m/s]
+EnableGravity = false
+
+[RANS]
+EddyViscosityModel = 3
+TurbulentSchmidtNumber = 0.7
+TurbulentPrandtlNumber = 0.85
+
+[Newton]
+MaxSteps = 10
+TargetSteps = 8
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = true