[freeflow][nonisothermal] Add enthalpy transported by diffusive fluxes

dumux/freeflow/compositional/staggered/fluxvariables.hh

@@ -54,14 +54,10 @@ class FreeflowNCFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
     using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
 
-    using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
-
+public:
     static constexpr auto numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
-
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-
-
-public:
+    using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
 
     /*!
     * \brief Computes the flux for the cell center residual.

dumux/freeflow/compositional/volumevariables.hh

@@ -217,6 +217,16 @@ public:
         return fluidState_.molarDensity(fluidSystemPhaseIdx);
     }
 
+    /*!
+     * \brief Returns the molar mass of a given component.
+     *
+     * \param compIdx the index of the component
+     */
+    Scalar molarMass(int compIdx) const
+    {
+        return FluidSystem::molarMass(compIdx);
+    }
+
      /*!
      * \brief Returns the diffusion coefficient \f$\mathrm{[m^2/s]}\f$
      *

dumux/freeflow/navierstokes/staggered/localresidual.hh

@@ -80,7 +80,7 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
 public:
-    using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance()>;
+    using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance(), (ModelTraits::numComponents() > 1)>;
 
     // account for the offset of the cell center privars within the PrimaryVariables container
     static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;

dumux/freeflow/nonisothermal/localresidual.hh

@@ -29,7 +29,7 @@
 namespace Dumux {
 
 // forward declaration
-template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod, bool enableEneryBalance>
+template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod, bool enableEneryBalance, bool isCompositional>
 class FreeFlowEnergyLocalResidualImplementation;
 
 /*!
@@ -37,19 +37,19 @@ class FreeFlowEnergyLocalResidualImplementation;
  * \brief Element-wise calculation of the local residual for non-isothermal
  *        free-flow models
  */
-template<class FVGridGeometry, class FluxVariables, bool enableEneryBalance>
+template<class FVGridGeometry, class FluxVariables, bool enableEneryBalance, bool isCompositional>
 using FreeFlowEnergyLocalResidual =
       FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
                                                 FluxVariables,
                                                 FVGridGeometry::discMethod,
-                                                enableEneryBalance>;
+                                                enableEneryBalance, isCompositional>;
 
 /*!
  * \ingroup FreeflowNIModel
  * \brief Specialization for isothermal models, does nothing
  */
-template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod>
-class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry, FluxVariables, discMethod, false>
+template<class FVGridGeometry, class FluxVariables, DiscretizationMethod discMethod, bool isCompositional>
+class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry, FluxVariables, discMethod, false, isCompositional>
 {
 public:
 
@@ -58,11 +58,6 @@ public:
     static void fluidPhaseStorage(Args&&... args)
     {}
 
-    //! do nothing for the isothermal case
-    template <typename... Args>
-    static void heatConvectionFlux(Args&&... args)
-    {}
-
     //! do nothing for the isothermal case
     template <typename... Args>
     static void heatFlux(Args&&... args)
@@ -71,18 +66,17 @@ public:
 
 /*!
  * \ingroup FreeflowNIModel
- * \brief Specialization for staggered non-isothermal models
+ * \brief Specialization for staggered one-phase, non-isothermal models
  */
 template<class FVGridGeometry, class FluxVariables>
 class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
                                                 FluxVariables,
                                                 DiscretizationMethod::staggered,
-                                                true>
+                                                true, false>
 {
     using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using HeatConductionType = typename FluxVariables::HeatConductionType;
 
 public:
 
@@ -123,10 +117,61 @@ public:
                                                                             upwindTerm,
                                                                             isOutflow);
         if(!isOutflow)
-            flux[localEnergyBalanceIdx] += HeatConductionType::flux(element,
-                                                               fvGeometry,
-                                                               elemVolVars,
-                                                               scvf);
+            flux[localEnergyBalanceIdx] += FluxVariables::HeatConductionType::flux(element,
+                                                                                   fvGeometry,
+                                                                                   elemVolVars,
+                                                                                   scvf);
+    }
+};
+
+/*!
+ * \ingroup FreeflowNIModel
+ * \brief Specialization for staggered compositional, non-isothermal models
+ */
+template<class FVGridGeometry, class FluxVariables>
+class FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
+                                                FluxVariables,
+                                                DiscretizationMethod::staggered,
+                                                true, true>
+    : public FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
+                                                       FluxVariables,
+                                                       DiscretizationMethod::staggered,
+                                                       true, false>
+{
+    using ParentType = FreeFlowEnergyLocalResidualImplementation<FVGridGeometry,
+                                                                 FluxVariables,
+                                                                 DiscretizationMethod::staggered,
+                                                                 true, false>;
+    using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+
+public:
+    //! The convective and conductive heat fluxes in the fluid phase
+    template<class NumEqVector, class Problem, class ElementVolumeVariables, class ElementFaceVariables>
+    static void heatFlux(NumEqVector& flux,
+                         const Problem& problem,
+                         const Element &element,
+                         const FVElementGeometry& fvGeometry,
+                         const ElementVolumeVariables& elemVolVars,
+                         const ElementFaceVariables& elemFaceVars,
+                         const SubControlVolumeFace& scvf)
+    {
+        ParentType::heatFlux(flux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+        static constexpr auto localEnergyBalanceIdx = NumEqVector::dimension - 1;
+        NumEqVector diffusiveFlux = FluxVariables::MolecularDiffusionType::flux(problem, fvGeometry, elemVolVars, scvf);
+        for (int compIdx = 0; compIdx < FluxVariables::numComponents; ++compIdx)
+        {
+            const bool insideIsUpstream = scvf.directionSign() == sign(diffusiveFlux[compIdx]);
+            const auto& upstreamVolVars = insideIsUpstream ? elemVolVars[scvf.insideScvIdx()] : elemVolVars[scvf.outsideScvIdx()];
+
+            // always use a mass-based calculation for the energy balance
+            if (FluxVariables::useMoles)
+                diffusiveFlux[compIdx] *= elemVolVars[scvf.insideScvIdx()].molarMass(compIdx);
+
+            flux[localEnergyBalanceIdx] += diffusiveFlux[compIdx] * upstreamVolVars.componentEnthalpy(compIdx);
+        }
     }
 };
 

dumux/freeflow/volumevariables.hh

@@ -179,8 +179,6 @@ public:
     /*!
      * \brief Returns the total internal energy of a phase in the
      *        sub-control volume.
-     *
-     * \param fluidSystemPhaseIdx The phase index
      */
     Scalar internalEnergy() const
     { return ParentType::asImp_().fluidState().internalEnergy(fluidSystemPhaseIdx); }
@@ -188,12 +186,16 @@ public:
     /*!
      * \brief Returns the total enthalpy of a phase in the sub-control
      *        volume.
-     *
-     * \param fluidSystemPhaseIdx The phase index
      */
     Scalar enthalpy() const
     { return ParentType::asImp_().fluidState().enthalpy(fluidSystemPhaseIdx); }
 
+    /*!
+     * \brief Returns the component enthalpy \f$\mathrm{[J/(kg*K)]}\f$ in the sub-control volume.
+     */
+    Scalar componentEnthalpy(unsigned int compIdx) const
+    { return FluidSystem::componentEnthalpy(ParentType::asImp_().fluidState(), fluidSystemPhaseIdx, compIdx); }
+
     /*!
      * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
      *        of the fluid phase in the sub-control volume.

dumux/material/fluidsystems/h2oair.hh

@@ -686,7 +686,8 @@ public:
 
         if (phaseIdx == phase0Idx)
         {
-            DUNE_THROW(Dune::NotImplemented, "The component enthalpies in the liquid phase are not implemented.");
+            // the liquid enthalpy is constant
+            return H2O::liquidEnthalpy(T, p);
         }
         else if (phaseIdx == phase1Idx)
         {

test/freeflow/ransnc/CMakeLists.txt

@@ -28,6 +28,6 @@ dune_add_test(NAME test_channel_lowrekepsilon2cni
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_channel_lowrekepsilon2cni.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_channel2cni-00033.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_channel2cni-00041.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_lowrekepsilon2cni test_channel2cni.input")
 target_compile_definitions(test_channel_lowrekepsilon2cni PUBLIC "LOWREKEPSILON=1" "NONISOTHERMAL=1")