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

---
 .../lowrekepsilon/staggered/fluxvariables.hh  | 20 +++++++++----------
 .../lowrekepsilon/staggered/localresidual.hh  | 16 +++++++--------
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
index a080d51cac..4ce64ab70d 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/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.dissipationTilde();
+            return volVars.dissipationTilde() * volVars.density();
         };
 
         flux[turbulentKineticEnergyEqIdx]
@@ -116,14 +116,14 @@ public:
         const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
         // effective diffusion coefficients
-        Scalar insideCoeff_k = insideVolVars.kinematicViscosity()
-                               + insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaK();
-        Scalar outsideCoeff_k = outsideVolVars.kinematicViscosity()
-                                + outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaK();
-        Scalar insideCoeff_e = insideVolVars.kinematicViscosity()
-                               + insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaEpsilon();
-        Scalar outsideCoeff_e = outsideVolVars.kinematicViscosity()
-                                + outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaEpsilon();
+        Scalar insideCoeff_k = insideVolVars.viscosity() + insideVolVars.kinematicEddyViscosity()
+                                                         * insideVolVars.density() / insideVolVars.sigmaK();
+        Scalar outsideCoeff_k = outsideVolVars.viscosity() + outsideVolVars.kinematicEddyViscosity()
+                                                           * outsideVolVars.density() / outsideVolVars.sigmaK();
+        Scalar insideCoeff_e = insideVolVars.viscosity() + insideVolVars.kinematicEddyViscosity()
+                                                         * insideVolVars.density() / insideVolVars.sigmaEpsilon();
+        Scalar outsideCoeff_e = outsideVolVars.viscosity() + outsideVolVars.kinematicEddyViscosity()
+                                                           * outsideVolVars.density() / outsideVolVars.sigmaEpsilon();
 
         // scale by extrusion factor
         insideCoeff_k *= insideVolVars.extrusionFactor();
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh
index b29e5bec58..1c692aa425 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh
@@ -82,8 +82,8 @@ public:
     {
         CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
 
-        storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
-        storage[dissipationEqIdx] = volVars.dissipationTilde();
+        storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy() * volVars.density();
+        storage[dissipationEqIdx] = volVars.dissipationTilde() * volVars.density();
 
         return storage;
     }
@@ -101,22 +101,22 @@ public:
         const auto& volVars = elemVolVars[scv];
 
         // production
-        source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.kinematicEddyViscosity()
+        source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.dynamicEddyViscosity()
                                                * volVars.stressTensorScalarProduct();
         source[dissipationEqIdx] += volVars.cOneEpsilon() * volVars.fOne()
                                     * volVars.dissipationTilde() / volVars.turbulentKineticEnergy()
-                                    * 2.0 * volVars.kinematicEddyViscosity()
+                                    * 2.0 * volVars.dynamicEddyViscosity()
                                     * volVars.stressTensorScalarProduct();
 
         // destruction
-        source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde();
-        source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo()
+        source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde() * volVars.density();
+        source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo() * volVars.density()
                                     * volVars.dissipationTilde() * volVars.dissipationTilde()
                                     / volVars.turbulentKineticEnergy();
 
         // dampening functions
-        source[turbulentKineticEnergyEqIdx] -= volVars.dValue();
-        source[dissipationEqIdx] += volVars.eValue();
+        source[turbulentKineticEnergyEqIdx] -= volVars.dValue() * volVars.density();
+        source[dissipationEqIdx] += volVars.eValue() * volVars.density();
 
         return source;
     }
-- 
GitLab