[freeflow][kepsilon] Implement wall functions for Scalar quantities

doc/handbook/dumux-handbook.bib

@@ -1714,3 +1714,14 @@ url={http://dx.doi.org/10.1007/s11242-015-0599-1}
   doi =       {10.2514/1.36541},
   url =       {https://doi.org/10.2514/1.36541}
 }
+
+@Book{Versteeg2009a,
+  title =     {{An Introduction to Computational Fluid Dynamics}},
+  publisher = {Pearson Education},
+  year =      {2009},
+  author =    {Versteeg, Henk and Malalasekra, Weeratunge},
+  isbn =      {978-0-13-127498-3},
+  pages =     {XII, 503},
+  address =   {Harlow},
+  edition =   {2},
+}

dumux/freeflow/navierstokes/problem.hh

@@ -76,6 +76,7 @@ class NavierStokesProblem : public NavierStokesParentProblem<TypeTag>
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
     using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
 
     enum {
@@ -190,6 +191,14 @@ public:
                                       const SubControlVolumeFace& localSubFace) const
     { return FacePrimaryVariables(0.0); }
 
+    //! \brief Returns an additional wall function flux for cell-centered quantities (only needed for RANS models)
+    CellCenterPrimaryVariables wallFunction(const Element& element,
+                                            const FVElementGeometry& fvGeometry,
+                                            const ElementVolumeVariables& elemVolVars,
+                                            const ElementFaceVariables& elemFaceVars,
+                                            const SubControlVolumeFace& scvf) const
+    { return CellCenterPrimaryVariables(0.0); }
+
 private:
 
     //! Returns the implementation of the problem (i.e. static polymorphism)

dumux/freeflow/navierstokes/staggered/localresidual.hh

@@ -79,6 +79,9 @@ class NavierStokesResidualImpl<TypeTag, DiscretizationMethod::staggered>
 
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
+    static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();
+    static constexpr bool isCompositional = ModelTraits::numComponents() > 1;
+
 public:
     using EnergyLocalResidual = FreeFlowEnergyLocalResidual<FVGridGeometry, FluxVariables, ModelTraits::enableEnergyBalance(), (ModelTraits::numComponents() > 1)>;
 
@@ -222,6 +225,7 @@ protected:
             if (scvf.boundary())
             {
                 const auto bcTypes = problem.boundaryTypes(element, scvf);
+                const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
 
                 // treat Dirichlet and outflow BCs
                 FluxVariables fluxVars;
@@ -231,19 +235,27 @@ protected:
                 EnergyLocalResidual::heatFlux(boundaryFlux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf);
 
                 // treat Neumann BCs, i.e. overwrite certain fluxes by user-specified values
+                static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
                 if(bcTypes.hasNeumann())
                 {
-                    static constexpr auto numEqCellCenter = CellCenterResidual::dimension;
                     for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
                     {
                         if(bcTypes.isNeumann(eqIdx + cellCenterOffset))
                         {
-                            const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
                             boundaryFlux[eqIdx] = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx + cellCenterOffset]
-                                                  * extrusionFactor * scvf.area();
+                                                                  * extrusionFactor * scvf.area();
                         }
                     }
                 }
+                for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+                {
+                    // use a wall function
+                    if(problem.useWallFunction(element, scvf, eqIdx + cellCenterOffset))
+                    {
+                        boundaryFlux[eqIdx] = problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[eqIdx]
+                                                                   * extrusionFactor * scvf.area();
+                    }
+                }
 
                 residual += boundaryFlux;
 

dumux/freeflow/rans/twoeq/kepsilon/problem.hh

@@ -68,9 +68,20 @@ class KEpsilonProblem : public RANSProblem<TypeTag>
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
     using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+    using Indices = typename ModelTraits::Indices;
+
+    static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();
+    static constexpr bool isCompositional = ModelTraits::numComponents() > 1;
+
+    // account for the offset of the cell center privars within the PrimaryVariables container
+    static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension;
+    static_assert(cellCenterOffset == ModelTraits::dim(), "cellCenterOffset must equal dim for staggered NavierStokes");
 
 public:
+    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
     //! The constructor sets the gravity, if desired by the user.
     KEpsilonProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
     : ParentType(fvGridGeometry)
@@ -319,6 +330,119 @@ public:
                                     * elemVolVars[scvf.insideScvIdx()].density());
     }
 
+    //! \brief Returns the flux for non-isothermal and compositional RANS models
+    template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+             typename std::enable_if_t<eB && compositional, int> = 0>
+    CellCenterPrimaryVariables wallFunction(const Element& element,
+                                            const FVElementGeometry& fvGeometry,
+                                            const ElementVolumeVariables& elemVolVars,
+                                            const ElementFaceVariables& elemFaceVars,
+                                            const SubControlVolumeFace& scvf) const
+    {
+        return wallFunctionComponent(element, fvGeometry, elemVolVars, elemFaceVars, scvf)
+               + wallFunctionEnergy(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+    }
+
+    //! \brief Returns the flux for isothermal and compositional RANS models
+    template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+             typename std::enable_if_t<!eB && compositional, int> = 0>
+    CellCenterPrimaryVariables wallFunction(const Element& element,
+                                            const FVElementGeometry& fvGeometry,
+                                            const ElementVolumeVariables& elemVolVars,
+                                            const ElementFaceVariables& elemFaceVars,
+                                            const SubControlVolumeFace& scvf) const
+    { return wallFunctionComponent(element, fvGeometry, elemVolVars, elemFaceVars, scvf); }
+
+    //! \brief Returns the flux for non-isothermal RANS models
+    template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+             typename std::enable_if_t<eB && !compositional, int> = 0>
+    CellCenterPrimaryVariables wallFunction(const Element& element,
+                                            const FVElementGeometry& fvGeometry,
+                                            const ElementVolumeVariables& elemVolVars,
+                                            const ElementFaceVariables& elemFaceVars,
+                                            const SubControlVolumeFace& scvf) const
+    { return wallFunctionEnergy(element, fvGeometry, elemVolVars, elemFaceVars, scvf); }
+
+    //! \brief Returns the flux for isothermal RANS models
+    template<bool eB = enableEnergyBalance, bool compositional = isCompositional,
+             typename std::enable_if_t<!eB && !compositional, int> = 0>
+    CellCenterPrimaryVariables wallFunction(const Element& element,
+                                            const FVElementGeometry& fvGeometry,
+                                            const ElementVolumeVariables& elemVolVars,
+                                            const ElementFaceVariables& elemFaceVars,
+                                            const SubControlVolumeFace& scvf) const
+    { return CellCenterPrimaryVariables(0.0); }
+
+    //! \brief Returns the component wall-function flux
+    CellCenterPrimaryVariables wallFunctionComponent(const Element& element,
+                                                     const FVElementGeometry& fvGeometry,
+                                                     const ElementVolumeVariables& elemVolVars,
+                                                     const ElementFaceVariables& elemFaceVars,
+                                                     const SubControlVolumeFace& scvf) const
+    {
+        auto wallFunctionFlux = CellCenterPrimaryVariables(0.0);
+        unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element);
+
+        // component mass fluxes
+        for (int compIdx = 0; compIdx < ModelTraits::numComponents(); ++compIdx)
+        {
+            if (Indices::replaceCompEqIdx == compIdx)
+                continue;
+
+            Scalar schmidtNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity()
+                                   / elemVolVars[scvf.insideScvIdx()].diffusionCoefficient(compIdx);
+            Scalar massConversionFactor = useMoles ? 1.0
+                                                   : FluidSystem::molarMass(compIdx);
+            wallFunctionFlux[compIdx] +=
+                -1.0 * (asImp_().dirichlet(element, scvf)[Indices::conti0EqIdx + compIdx]
+                        - elemVolVars[scvf.insideScvIdx()].moleFraction(compIdx))
+                * elemVolVars[scvf.insideScvIdx()].molarDensity()
+                * uStarNominal(elementID)
+                / asImp_().turbulentSchmidtNumber()
+                / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementID) * 9.793)
+                    + pFunction(schmidtNumber, asImp_().turbulentSchmidtNumber()));
+        }
+
+        return wallFunctionFlux;
+    }
+
+    //! \brief Returns the energy wall-function flux
+    CellCenterPrimaryVariables wallFunctionEnergy(const Element& element,
+                                                  const FVElementGeometry& fvGeometry,
+                                                  const ElementVolumeVariables& elemVolVars,
+                                                  const ElementFaceVariables& elemFaceVars,
+                                                  const SubControlVolumeFace& scvf) const
+    {
+        auto wallFunctionFlux = CellCenterPrimaryVariables(0.0);
+        unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element);
+        // energy fluxes
+        Scalar prandtlNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity()
+                               * elemVolVars[scvf.insideScvIdx()].density()
+                               * elemVolVars[scvf.insideScvIdx()].heatCapacity()
+                               / elemVolVars[scvf.insideScvIdx()].thermalConductivity();
+        wallFunctionFlux[Indices::energyBalanceIdx - cellCenterOffset] +=
+            -1.0 * (asImp_().dirichlet(element, scvf)[Indices::temperatureIdx]
+                    - elemVolVars[scvf.insideScvIdx()].temperature())
+            * elemVolVars[scvf.insideScvIdx()].density()
+            * elemVolVars[scvf.insideScvIdx()].heatCapacity()
+            * uStarNominal(elementID)
+            / asImp_().turbulentPrandtlNumber()
+            / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementID) * 9.793)
+                + pFunction(prandtlNumber, asImp_().turbulentPrandtlNumber()));
+
+        return wallFunctionFlux;
+    }
+
+    //! \brief Returns the value of the P-function after Jayatilleke \cite Versteeg2009a
+    const Scalar pFunction(Scalar molecularNumber, Scalar turbulentNumber) const
+    {
+        using std::pow;
+        using std::exp;
+        return 9.24
+               * (pow(molecularNumber / turbulentNumber, 0.75) - 1.0)
+               * (1.0 + 0.28 * exp(-0.007 * molecularNumber / turbulentNumber));
+    }
+
     //! \brief Returns the \$f C_{\mu} \$f constant
     const Scalar cMu() const
     {