[ransnc][oneeq] Implement a compositional one-equation

dumux/freeflow/compositional/oneeqncmodel.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 FreeflowNCModel
+ *
+ * \brief A single-phase, multi-component one-equation model
+ *
+ * \copydoc Dumux::FreeflowNCModel
+ */
+
+#ifndef DUMUX_ONEEQ_NC_MODEL_HH
+#define DUMUX_ONEEQ_NC_MODEL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+#include <dumux/freeflow/nonisothermal/vtkoutputfields.hh>
+#include <dumux/freeflow/rans/oneeq/model.hh>
+
+#include "vtkoutputfields.hh"
+
+namespace Dumux {
+
+///////////////////////////////////////////////////////////////////////////
+// properties for the single-phase, multi-component one-equation model
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the single-phase, multi-component isothermal one-equation model
+NEW_TYPE_TAG(OneEqNC, INHERITS_FROM(NavierStokesNC));
+
+///////////////////////////////////////////////////////////////////////////
+// default property values
+///////////////////////////////////////////////////////////////////////////
+
+/*!
+ * \ingroup FreeflowNCModel
+ * \brief Traits for the one-equation multi-component model
+ */
+template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
+struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+{
+    //! There are as many momentum balance equations as dimensions
+    //! and as many balance equations as components.
+    static constexpr int numEq() { return dimension+nComp+1; }
+
+    //! The model does include a turbulence model
+    static constexpr bool usesTurbulenceModel() { return true; }
+
+    //! the indices
+    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, OneEqIndices<dimension, nComp>>;
+};
+
+//!< states some specifics of the isothermal multi-component one-equation model
+SET_PROP(OneEqNC, ModelTraits)
+{
+private:
+    using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
+    static constexpr int dimension = GridView::dimension;
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static constexpr int numComponents = FluidSystem::numComponents;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+public:
+    using type = OneEqNCModelTraits<dimension, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+};
+
+//! Set the volume variables property
+SET_PROP(OneEqNC, 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 NCVolVars = FreeflowNCVolumeVariables<Traits>;
+public:
+    using type = OneEqVolumeVariables<Traits, NCVolVars>;
+};
+
+//! The local residual
+SET_PROP(OneEqNC, LocalResidual)
+{
+private:
+    using BaseLocalResidual = FreeflowNCResidual<TypeTag>;
+public:
+    using type = OneEqResidual<TypeTag, BaseLocalResidual>;
+};
+
+//! The flux variables
+SET_PROP(OneEqNC, FluxVariables)
+{
+private:
+    using BaseFluxVariables = FreeflowNCFluxVariables<TypeTag>;
+public:
+    using type = OneEqFluxVariables<TypeTag, BaseFluxVariables>;
+};
+
+//! The specific vtk output fields
+SET_PROP(OneEqNC, VtkOutputFields)
+{
+private:
+    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+    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 SinglePhaseVtkOutputFields = OneEqVtkOutputFields<FVGridGeometry>;
+public:
+    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+};
+
+//////////////////////////////////////////////////////////////////////////
+// Property values for non-isothermal multi-component one-equation model
+//////////////////////////////////////////////////////////////////////////
+
+//! The type tags for the single-phase, multi-component non-isothermal one-equation models
+NEW_TYPE_TAG(OneEqNCNI, INHERITS_FROM(NavierStokesNCNI));
+
+//! The model traits of the non-isothermal model
+SET_PROP(OneEqNCNI, 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;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    using IsothermalModelTraits = OneEqNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+public:
+    using type = FreeflowNIModelTraits<IsothermalModelTraits>;
+};
+
+//! Set the volume variables property
+SET_PROP(OneEqNCNI, 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 NCVolVars = FreeflowNCVolumeVariables<Traits>;
+public:
+    using type = OneEqVolumeVariables<Traits, NCVolVars>;
+};
+
+//! The local residual
+SET_PROP(OneEqNCNI, LocalResidual)
+{
+private:
+    using BaseLocalResidual = FreeflowNCResidual<TypeTag>;
+public:
+    using type = OneEqResidual<TypeTag, BaseLocalResidual>;
+};
+
+//! The flux variables
+SET_PROP(OneEqNCNI, FluxVariables)
+{
+private:
+    using BaseFluxVariables = FreeflowNCFluxVariables<TypeTag>;
+public:
+    using type = OneEqFluxVariables<TypeTag, BaseFluxVariables>;
+};
+
+//! The specific vtk output fields
+SET_PROP(OneEqNCNI, VtkOutputFields)
+{
+private:
+    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+    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 BaseVtkOutputFields = OneEqVtkOutputFields<FVGridGeometry>;
+    using NonIsothermalFields = FreeflowNonIsothermalVtkOutputFields<BaseVtkOutputFields, ModelTraits>;
+public:
+    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+};
+
+// \}
+} // end namespace Properties
+} // end namespace Dumux
+
+#endif

test/freeflow/ransnc/CMakeLists.txt

@@ -60,6 +60,21 @@ dune_add_test(NAME test_flatplate_lowrekepsilon2cni
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_flatplate_lowrekepsilon2cni test_flatplate2cni.input")
 target_compile_definitions(test_flatplate_lowrekepsilon2cni PUBLIC "LOWREKEPSILON=1" "NONISOTHERMAL=1")
 
+dune_add_test(NAME test_flatplate_oneeq2c
+              SOURCES test_flatplate.cc
+              COMPILE_ONLY CMAKE_GUARD HAVE_UMFPACK)
+target_compile_definitions(test_flatplate_oneeq2c PUBLIC "ONEEQ=1")
+
+dune_add_test(NAME test_flatplate_oneeq2cni
+              SOURCES test_flatplate.cc
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS       --script fuzzy
+                             --files ${CMAKE_SOURCE_DIR}/test/references/test_flatplate_oneeq2cni.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_flatplate2cni-00034.vtu
+                             --command "${CMAKE_CURRENT_BINARY_DIR}/test_flatplate_oneeq2cni test_flatplate2cni.input")
+target_compile_definitions(test_flatplate_oneeq2cni PUBLIC "ONEEQ=1" "NONISOTHERMAL=1")
+
 dune_add_test(NAME test_flatplate_zeroeq2c
               SOURCES test_flatplate.cc
               CMAKE_GUARD HAVE_UMFPACK

test/freeflow/ransnc/flatplatetestproblem.hh

@@ -39,6 +39,9 @@
 #elif KOMEGA
 #include <dumux/freeflow/compositional/komegancmodel.hh>
 #include <dumux/freeflow/rans/twoeq/komega/problem.hh>
+#elif ONEEQ
+#include <dumux/freeflow/compositional/oneeqncmodel.hh>
+#include <dumux/freeflow/rans/oneeq/problem.hh>
 #else
 #include <dumux/freeflow/compositional/zeroeqncmodel.hh>
 #include <dumux/freeflow/rans/zeroeq/problem.hh>
@@ -59,6 +62,8 @@ namespace Properties
   NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KEpsilonNCNI));
   #elif KOMEGA
   NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KOmegaNCNI));
+  #elif ONEEQ
+  NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNCNI));
   #else
   NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNCNI));
   #endif
@@ -69,6 +74,8 @@ namespace Properties
   NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KEpsilonNC));
   #elif KOMEGA
   NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, KOmegaNC));
+  #elif ONEEQ
+  NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, OneEqNC));
   #else
   NEW_TYPE_TAG(FlatPlateNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNC));
   #endif
@@ -123,6 +130,10 @@ class FlatPlateNCTestProblem : public KEpsilonProblem<TypeTag>
 class FlatPlateNCTestProblem : public KOmegaProblem<TypeTag>
 {
     using ParentType = KOmegaProblem<TypeTag>;
+#elif ONEEQ
+class FlatPlateNCTestProblem : public OneEqProblem<TypeTag>
+{
+    using ParentType = OneEqProblem<TypeTag>;
 #else
 class FlatPlateNCTestProblem : public ZeroEqProblem<TypeTag>
 {
@@ -166,6 +177,7 @@ public:
         Scalar density = FluidSystem::density(fluidState, phaseIdx);
         Scalar kinematicViscosity = FluidSystem::viscosity(fluidState, phaseIdx) / density;
         Scalar diameter = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
+        viscosityTilde_ = 1e-3 * turbulenceProperties.viscosityTilde(inletVelocity_, diameter, kinematicViscosity);
         turbulentKineticEnergy_ = turbulenceProperties.turbulentKineticEnergy(inletVelocity_, diameter, kinematicViscosity);
 #if KOMEGA
         dissipation_ = turbulenceProperties.dissipationRate(inletVelocity_, diameter, kinematicViscosity);
@@ -231,6 +243,9 @@ public:
 #if KEPSILON || KOMEGA || LOWREKEPSILON
             values.setDirichlet(Indices::turbulentKineticEnergyIdx);
             values.setDirichlet(Indices::dissipationIdx);
+#endif
+#if ONEEQ
+            values.setDirichlet(Indices::viscosityTildeIdx);
 #endif
         }
         else if(isOutlet_(globalPos))
@@ -245,6 +260,9 @@ public:
 #if KEPSILON || KOMEGA || LOWREKEPSILON
             values.setOutflow(Indices::turbulentKineticEnergyEqIdx);
             values.setOutflow(Indices::dissipationEqIdx);
+#endif
+#if ONEEQ
+            values.setOutflow(Indices::viscosityTildeIdx);
 #endif
         }
         else if(isOnWall(globalPos))
@@ -263,6 +281,9 @@ public:
 #elif KOMEGA
             values.setDirichlet(Indices::turbulentKineticEnergyEqIdx);
             values.setDirichletCell(Indices::dissipationIdx);
+#endif
+#if ONEEQ
+            values.setDirichlet(Indices::viscosityTildeIdx);
 #endif
         }
         else
@@ -362,6 +383,14 @@ public:
         }
 #endif
 
+#if ONEEQ
+        values[Indices::viscosityTildeIdx] = viscosityTilde_;
+        if (isOnWall(globalPos))
+        {
+            values[Indices::viscosityTildeIdx] = 0.0;
+        }
+#endif
+
         return values;
     }
 
@@ -397,6 +426,7 @@ private:
 
     const Scalar eps_;
     Scalar inletVelocity_;
+    Scalar viscosityTilde_;
     Scalar turbulentKineticEnergy_;
     Scalar dissipation_;
     TimeLoopPtr timeLoop_;