[freeflow] change diffusive flux computation to mass averaged reference

system and adapt units of fluxes accordingly

dumux/freeflow/compositional/navierstokesncmodel.hh

@@ -26,7 +26,7 @@
  * \f[
  *    \frac{\partial \left(\varrho X^\kappa\right)}{\partial t}
  *    + \nabla \cdot \left( \varrho {\boldsymbol{v}} X^\kappa
- *    - (D^\kappa + D_\text{t}) \varrho \frac{M^\kappa}{M} \textbf{grad}\, x^\kappa \right)
+ *    - (D^\kappa + D_\text{t}) \varrho \textbf{grad}\, X^\kappa \right)
  *    - q^\kappa = 0
  * \f]
  *
@@ -35,7 +35,7 @@
  *    \frac{\partial \varrho_g}{\partial t}
  *    + \nabla \cdot \left(
  *        \varrho {\boldsymbol{v}}
- *        - \sum_\kappa (D^\kappa + D_\text{t}) \varrho \frac{M^\kappa}{M} \textbf{grad}\, x^\kappa
+ *        - \sum_\kappa (D^\kappa + D_\text{t}) \varrho \textbf{grad}\, X^\kappa
  *      \right)
  *    - q = 0
  * \f]

dumux/freeflow/compositional/staggered/fluxvariables.hh

@@ -86,7 +86,13 @@ public:
 
             flux[compIdx] = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTerm);
 
-            flux[compIdx] += useMoles ? diffusiveFluxes[compIdx] : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+            //check for the reference system and adapt units of the diffusive flux accordingly.
+            if (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+            {
+                flux[compIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
+            }
+            else
+                flux[compIdx] += useMoles ? diffusiveFluxes[compIdx] : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
         }
 
         // in case one balance is substituted by the total mass balance

dumux/freeflow/compositional/volumevariables.hh

@@ -245,6 +245,14 @@ public:
         return FluidSystem::molarMass(compIdx);
     }
 
+    /*!
+     * \brief Returns the average molar mass \f$\mathrm{[kg/mol]}\f$$ the of the fluid phase.
+     *
+     * \param phaseIdx The phase index
+     */
+    Scalar averageMolarMass(const int phaseIdx = 0) const
+    { return fluidState_.averageMolarMass(phaseIdx); }
+
    /*!
     * \brief Returns the diffusion coefficient \f$\mathrm{[m^2/s]}\f$
     *

dumux/freeflow/nonisothermal/localresidual.hh

@@ -25,6 +25,7 @@
 #define DUMUX_FREE_FLOW_ENERGY_LOCAL_RESIDUAL_HH
 
 #include <dumux/discretization/method.hh>
+#include <dumux/flux/referencesystem.hh>
 
 namespace Dumux {
 
@@ -158,11 +159,11 @@ public:
             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);
+            if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+                flux[localEnergyBalanceIdx] += diffusiveFlux[compIdx] * upstreamVolVars.componentEnthalpy(compIdx);
+            else
+                flux[localEnergyBalanceIdx] += diffusiveFlux[compIdx] * upstreamVolVars.componentEnthalpy(compIdx)* elemVolVars[scvf.insideScvIdx()].molarMass(compIdx);
 
-            flux[localEnergyBalanceIdx] += diffusiveFlux[compIdx] * upstreamVolVars.componentEnthalpy(compIdx);
         }
     }
 };