From af178649b587214430343343d90d97659bf914cf Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Thu, 5 Sep 2019 15:22:44 +0200
Subject: [PATCH] [freeflow] change diffusive flux computation to mass averaged
 reference system and adapt units of fluxes accordingly

---
 dumux/freeflow/compositional/navierstokesncmodel.hh     | 4 ++--
 dumux/freeflow/compositional/staggered/fluxvariables.hh | 8 +++++++-
 dumux/freeflow/compositional/volumevariables.hh         | 8 ++++++++
 dumux/freeflow/nonisothermal/localresidual.hh           | 9 +++++----
 4 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/dumux/freeflow/compositional/navierstokesncmodel.hh b/dumux/freeflow/compositional/navierstokesncmodel.hh
index c909318257..b7958d0155 100644
--- a/dumux/freeflow/compositional/navierstokesncmodel.hh
+++ b/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]
diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh
index 32ae6f4610..de283bab6f 100644
--- a/dumux/freeflow/compositional/staggered/fluxvariables.hh
+++ b/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
diff --git a/dumux/freeflow/compositional/volumevariables.hh b/dumux/freeflow/compositional/volumevariables.hh
index a9ccb5ae87..8be88d881c 100644
--- a/dumux/freeflow/compositional/volumevariables.hh
+++ b/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$
     *
diff --git a/dumux/freeflow/nonisothermal/localresidual.hh b/dumux/freeflow/nonisothermal/localresidual.hh
index 773157340b..0a5475a3a7 100644
--- a/dumux/freeflow/nonisothermal/localresidual.hh
+++ b/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);
         }
     }
 };
-- 
GitLab