[rans][lowrekepsilon] add density in all terms of the lowrekepsilon model

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();

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;
     }