From 2a96db9240998020dd7ca698f1c15933cc56c6f3 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 6 Aug 2020 17:35:23 +0200
Subject: [PATCH] [rans][kepsilon] add density in all terms of the kepsilon
 model

---
 dumux/freeflow/rans/twoeq/kepsilon/problem.hh |  2 +-
 .../twoeq/kepsilon/staggered/fluxvariables.hh | 21 ++++++-------------
 .../twoeq/kepsilon/staggered/localresidual.hh | 12 +++++------
 3 files changed, 13 insertions(+), 22 deletions(-)

diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 3f73d9b800..5541f7aa23 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -325,7 +325,7 @@ public:
     {
         unsigned int elementIdx = asImp_().gridGeometry().elementMapper().index(element);
         return FacePrimaryVariables(asImp_().tangentialMomentumWallFunction(elementIdx, elemFaceVars[scvf].velocitySelf())
-                                    * elemVolVars[scvf.insideScvIdx()].density());
+                                    * asImp_().storedDensity(elementIdx) );
     }
 
     //! \brief Returns the flux for non-isothermal and compositional RANS models
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
index 5050c4ad10..0eca6c4ecc 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
@@ -97,11 +97,11 @@ public:
         // calculate advective flux
         auto upwindTermK = [](const auto& volVars)
         {
-            return volVars.turbulentKineticEnergy();
+            return volVars.turbulentKineticEnergy() * volVars.density();
         };
         auto upwindTermEpsilon = [](const auto& volVars)
         {
-            return volVars.dissipation();
+            return volVars.dissipation() * volVars.density();
         };
 
         flux[turbulentKineticEnergyEqIdx]
@@ -116,19 +116,10 @@ public:
         const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
         // effective diffusion coefficients
-        Scalar insideCoeff_k = insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaK();
-        Scalar outsideCoeff_k = outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaK();
-        Scalar insideCoeff_e = insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaEpsilon();
-        Scalar outsideCoeff_e = outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaEpsilon();
-        static const auto kEpsilonEnableKinematicViscosity_
-            = getParamFromGroup<bool>(problem.paramGroup(), "KEpsilon.EnableKinematicViscosity", true);
-        if (kEpsilonEnableKinematicViscosity_)
-        {
-            insideCoeff_k += insideVolVars.kinematicViscosity();
-            outsideCoeff_k += outsideVolVars.kinematicViscosity();
-            insideCoeff_e += insideVolVars.kinematicViscosity();
-            outsideCoeff_e += outsideVolVars.kinematicViscosity();
-        }
+        Scalar insideCoeff_k = (insideVolVars.dynamicEddyViscosity() / insideVolVars.sigmaK()) + insideVolVars.viscosity();
+        Scalar outsideCoeff_k = (outsideVolVars.dynamicEddyViscosity() / outsideVolVars.sigmaK()) + outsideVolVars.viscosity();
+        Scalar insideCoeff_e = (insideVolVars.dynamicEddyViscosity() / insideVolVars.sigmaEpsilon()) + insideVolVars.viscosity();
+        Scalar outsideCoeff_e = (outsideVolVars.dynamicEddyViscosity() / outsideVolVars.sigmaEpsilon()) + outsideVolVars.viscosity();
 
         // scale by extrusion factor
         insideCoeff_k *= insideVolVars.extrusionFactor();
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh
index acfa3e939d..df600f9226 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh
@@ -89,8 +89,8 @@ public:
     {
         CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
 
-        storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
-        storage[dissipationEqIdx] = volVars.dissipation();
+        storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy() * volVars.density();
+        storage[dissipationEqIdx] = volVars.dissipation() * volVars.density();
 
         return storage;
     }
@@ -113,23 +113,23 @@ public:
         // turbulence production is equal to dissipation -> exclude both terms (according to local equilibrium hypothesis, see FLUENT)
         if (!volVars.isMatchingPoint())
         {
-            source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.kinematicEddyViscosity()
+            source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.dynamicEddyViscosity()
                                                    * volVars.stressTensorScalarProduct();
         }
         source[dissipationEqIdx] += volVars.cOneEpsilon()
                                     * dissipation / turbulentKineticEnergy
-                                    * 2.0 * volVars.kinematicEddyViscosity()
+                                    * 2.0 * volVars.dynamicEddyViscosity()
                                     * volVars.stressTensorScalarProduct();
 
         // destruction
         // turbulence production is equal to dissipation -> exclude both terms (according to local equilibrium hypothesis, see FLUENT)
         if (!volVars.isMatchingPoint())
         {
-            source[turbulentKineticEnergyEqIdx] -= dissipation;
+            source[turbulentKineticEnergyEqIdx] -= dissipation * volVars.density();
         }
         source[dissipationEqIdx] -= volVars.cTwoEpsilon()
                                     * dissipation * dissipation
-                                    / turbulentKineticEnergy;
+                                    / turbulentKineticEnergy * volVars.density();
 
         return source;
     }
-- 
GitLab