[navierstokes] Allow use of 2p immiscible fluid systems for Stokes

* this makes coupling with, e.g., Darcy domains much easiert, as the same fluid system can be used for both domains
* add fluidSystemPhaseIdx to model traits and indices, where necessary
* add some docu

dumux/freeflow/compositional/kepsilonncmodel.hh

@@ -68,7 +68,7 @@ struct KEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phase
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, KEpsilonIndices<dimension, nComp>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, KEpsilonIndices<dimension, nComp, phaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model

dumux/freeflow/compositional/komegancmodel.hh

@@ -56,9 +56,15 @@ NEW_TYPE_TAG(KOmegaNC, INHERITS_FROM(NavierStokesNC));
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the low-Reynolds k-epsilon multi-component model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam nComp The number of components to be considered
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
+ * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
+ * \tparam useM Use molar or mass balances
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
-struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
+struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -68,7 +74,7 @@ struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseId
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, KOmegaIndices<dimension, nComp>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, KOmegaIndices<dimension, nComp, fluidSystemPhaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model

dumux/freeflow/compositional/lowrekepsilonncmodel.hh

@@ -56,9 +56,15 @@ NEW_TYPE_TAG(LowReKEpsilonNC, INHERITS_FROM(NavierStokesNC));
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the low-Reynolds k-epsilon multi-component model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam nComp The number of components to be considered
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
+ * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
+ * \tparam useM Use molar or mass balances
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
-struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
+struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -68,7 +74,7 @@ struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp,
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, LowReKEpsilonIndices<dimension, nComp>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, LowReKEpsilonIndices<dimension, nComp, fluidSystemPhaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model

dumux/freeflow/compositional/navierstokesncmodel.hh

@@ -76,9 +76,15 @@ namespace Dumux {
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the multi-component free-flow model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam nComp The number of components to be considered
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
+ * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
+ * \tparam useM Use molar or mass balances
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useM>
-struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useM>
+struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -94,7 +100,7 @@ struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension>
     static constexpr bool enableMolecularDiffusion() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, NavierStokesIndices<dimension>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, NavierStokesIndices<dimension, fluidSystemPhaseIdx>>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -127,6 +133,10 @@ private:
     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);
+
+    static_assert(phaseIdx >= 0 && phaseIdx < FluidSystem::numPhases,
+                  "PhaseIdx must be non-negative and smaller than the number of phases");
+
 public:
     using type = NavierStokesNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
 };

dumux/freeflow/compositional/oneeqncmodel.hh

@@ -56,9 +56,15 @@ NEW_TYPE_TAG(OneEqNC, INHERITS_FROM(NavierStokesNC));
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the one-equation multi-component model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam nComp The number of components to be considered
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
+ * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
+ * \tparam useM Use molar or mass balances
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
-struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
+struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -68,7 +74,7 @@ struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, OneEqIndices<dimension, nComp>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, OneEqIndices<dimension, nComp, fluidSystemPhaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component one-equation model

dumux/freeflow/navierstokes/indices.hh

@@ -31,8 +31,9 @@ namespace Dumux {
  * \brief The common indices for the isothermal Navier-Stokes model.
  *
  * \tparam dimension The dimension of the problem
+ * \tparam fsPhaseIdx The the index of the phase used for the fluid system
  */
-template <int dimension>
+template <int dimension, int fsPhaseIdx>
 struct NavierStokesIndices
 {
     static constexpr int dimXIdx = 0; //!< Index of the x-component of a vector of size dim
@@ -63,7 +64,7 @@ struct NavierStokesIndices
     }
 
     //! The index of the fluid phase in the fluid system (for compatibility reasons)
-    static constexpr int fluidSystemPhaseIdx = 0;
+    static constexpr int fluidSystemPhaseIdx = fsPhaseIdx;
 };
 
 } // end namespace Dumux

dumux/freeflow/navierstokes/model.hh

@@ -70,8 +70,11 @@ namespace Dumux {
 /*!
  * \ingroup NavierStokesModel
  * \brief Traits for the Navier-Stokes model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension>
+template<int dimension, int fluidSystemPhaseIdx>
 struct NavierStokesModelTraits
 {
     //! The dimension of the model
@@ -100,7 +103,7 @@ struct NavierStokesModelTraits
     static constexpr bool usesTurbulenceModel() { return false; }
 
     //! the indices
-    using Indices = NavierStokesIndices<dim()>;
+    using Indices = NavierStokesIndices<dim(), fluidSystemPhaseIdx>;
 };
 
 /*!
@@ -152,9 +155,13 @@ SET_PROP(NavierStokes, ModelTraits)
 {
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
-    static constexpr int dim = GridView::dimension;
+    static constexpr auto dim = GridView::dimension;
+    static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+    static_assert(phaseIdx >= 0 && phaseIdx < GET_PROP_TYPE(TypeTag, FluidSystem)::numPhases,
+                  "PhaseIdx must be non-negative and smaller than the number of phases");
 public:
-    using type = NavierStokesModelTraits<dim>;
+    using type = NavierStokesModelTraits<dim, phaseIdx>;
 };
 
 /*!
@@ -183,9 +190,8 @@ private:
     using FST = typename GET_PROP_TYPE(TypeTag, FluidState);
     using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
-    static_assert(MT::numPhases() == 1 && MT::numComponents() == 1 &&
-                  FSY::numPhases == 1 && FSY::numComponents == 1,
-                  "The Navier-Stokes model only works with a single-phase fluid system.");
+    static_assert(!FSY::isMiscible(),
+                  "The Navier-Stokes model only works with immiscible fluid systems.");
 
     using Traits = NavierStokesVolumeVariablesTraits<PV, FSY, FST, MT>;
 public:
@@ -215,8 +221,9 @@ SET_PROP(NavierStokesNI, ModelTraits)
 {
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
-    static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = NavierStokesModelTraits<dim>;
+    static constexpr auto dim = GridView::dimension;
+    static constexpr auto fluidSystemPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = NavierStokesModelTraits<dim, fluidSystemPhaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };

dumux/freeflow/rans/model.hh

@@ -65,9 +65,12 @@ SET_BOOL_PROP(RANS, EnableInertiaTerms, true); //!< Explicitly force the conside
 /*!
  * \ingroup RANSModel
  * \brief Traits for the Reynolds-averaged Navier-Stokes model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension>
-struct RANSModelTraits : NavierStokesModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct RANSModelTraits : NavierStokesModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The model does include a turbulence model
     static constexpr bool usesTurbulenceModel() { return true; }
@@ -79,8 +82,12 @@ SET_PROP(RANS, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+    static_assert(phaseIdx >= 0 && phaseIdx < GET_PROP_TYPE(TypeTag, FluidSystem)::numPhases,
+                  "PhaseIdx must be non-negative and smaller than the number of phases");
 public:
-    using type = RANSModelTraits<dim>;
+    using type = RANSModelTraits<dim, phaseIdx>;
 };
 
 //! The specific vtk output fields
@@ -105,7 +112,9 @@ SET_PROP(RANSNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = RANSModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+    using IsothermalTraits = RANSModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };

dumux/freeflow/rans/oneeq/indices.hh

@@ -35,9 +35,10 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents>
-struct OneEqIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct OneEqIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 public:
     static constexpr auto viscosityTildeEqIdx = dimension + numComponents;

dumux/freeflow/rans/oneeq/model.hh

@@ -94,9 +94,12 @@ namespace Properties {
 /*!
  * \ingroup OneEqModel
  * \brief Traits for the Spalart-Allmaras model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension>
-struct OneEqModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct OneEqModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -109,7 +112,7 @@ struct OneEqModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = OneEqIndices<dim(), numComponents()>;
+    using Indices = OneEqIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -125,8 +128,9 @@ SET_PROP(OneEq, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = OneEqModelTraits<dim>;
+    using type = OneEqModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -184,7 +188,8 @@ SET_PROP(OneEqNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = OneEqModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = OneEqModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };

dumux/freeflow/rans/twoeq/kepsilon/indices.hh

@@ -35,9 +35,10 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents>
-struct KEpsilonIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct KEpsilonIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 public:
     static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;

dumux/freeflow/rans/twoeq/kepsilon/model.hh

@@ -81,9 +81,12 @@ namespace Properties {
 /*!
  * \ingroup KEpsilonModel
  * \brief Traits for the k-epsilon model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension>
-struct KEpsilonModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct KEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -96,7 +99,7 @@ struct KEpsilonModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = KEpsilonIndices<dim(), numComponents()>;
+    using Indices = KEpsilonIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -112,8 +115,9 @@ SET_PROP(KEpsilon, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = KEpsilonModelTraits<dim>;
+    using type = KEpsilonModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -171,7 +175,8 @@ SET_PROP(KEpsilonNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = KEpsilonModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = KEpsilonModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };

dumux/freeflow/rans/twoeq/komega/indices.hh

@@ -35,9 +35,10 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents>
-struct KOmegaIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct KOmegaIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 public:
     static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;

dumux/freeflow/rans/twoeq/komega/model.hh

@@ -88,9 +88,12 @@ namespace Properties {
 /*!
  *\ingroup KOmegaModel
  * \brief Traits for the k-omega model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension>
-struct KOmegaModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct KOmegaModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -103,7 +106,7 @@ struct KOmegaModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! The indices
-    using Indices = KOmegaIndices<dim(), numComponents()>;
+    using Indices = KOmegaIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -119,8 +122,9 @@ SET_PROP(KOmega, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = KOmegaModelTraits<dim>;
+    using type = KOmegaModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -179,7 +183,8 @@ SET_PROP(KOmegaNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = KOmegaModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = KOmegaModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };

dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh

@@ -35,9 +35,10 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents>
-struct LowReKEpsilonIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct LowReKEpsilonIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 public:
     static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;

dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh

@@ -97,9 +97,12 @@ namespace Properties {
 /*!
  * \ingroup LowReKEpsilonModel
  * \brief Traits for the low-Reynolds k-epsilon model
+ *
+ * \tparam dimension The dimension of the problem
+ * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension>
-struct LowReKEpsilonModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct LowReKEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -112,7 +115,7 @@ struct LowReKEpsilonModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = LowReKEpsilonIndices<dim(), numComponents()>;
+    using Indices = LowReKEpsilonIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -128,8 +131,9 @@ SET_PROP(LowReKEpsilon, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = LowReKEpsilonModelTraits<dim>;
+    using type = LowReKEpsilonModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -187,7 +191,8 @@ SET_PROP(LowReKEpsilonNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = LowReKEpsilonModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = LowReKEpsilonModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };