[zeroeq] Implement first skeleton for more sophisticated turbulence model (modified vanDriest)

dumux/freeflow/rans/CMakeLists.txt

@@ -2,10 +2,6 @@ add_subdirectory("zeroeq")
 
 #install headers
 install(FILES
-# fluxvariables.hh
-# fluxvariablescache.hh
-# indices.hh
-# localresidual.hh
 model.hh
 problem.hh
 volumevariables.hh

dumux/freeflow/rans/problem.hh

@@ -66,6 +66,7 @@ class RANSProblem : public NavierStokesParentProblem<TypeTag>
       };
     // TODO: dim or dimWorld appropriate here?
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+    using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 
 public:
     //! The constructor sets the gravity, if desired by the user.
@@ -84,12 +85,18 @@ public:
      */
     void updateStaticWallProperties() const
     {
+        std::cout << "Update static wall properties. ";
+
         // update size and initial values of the global vectors
         wallElementIDs_.resize(this->fvGridGeometry().elementMapper().size());
         wallDistances_.resize(this->fvGridGeometry().elementMapper().size());
-            for (unsigned int i = 0; i < wallElementIDs_.size(); ++i)
+        velocityGradients_.resize(this->fvGridGeometry().elementMapper().size());
+        kinematicViscosity_.resize(this->fvGridGeometry().elementMapper().size());
+        for (unsigned int i = 0; i < wallElementIDs_.size(); ++i)
         {
             wallDistances_[i] = std::numeric_limits<Scalar>::max();
+            velocityGradients_[i] = DimMatrix(0.0);
+            kinematicViscosity_[i] = 0.0;
         }
 
         // retrieve all wall intersections and corresponding elements
@@ -108,7 +115,7 @@ public:
                 }
             }
         }
-        std::cout << "numWallIntersections: " << wallPositions.size() << std::endl;
+        std::cout << "NumWallIntersections=" << wallPositions.size() << std::endl;
 
         // search for shortest distance to wall for each element
         for (const auto& element : elements(gridView))
@@ -132,7 +139,17 @@ public:
      */
     void updateDynamicWallProperties() const
     {
-        return;
+        std::cout << "Update dynamic wall properties." << std::endl;
+
+        auto& gridView(this->fvGridGeometry().gridView());
+        for (const auto& element : elements(gridView))
+        {
+            // TODO call calculate velocity gradients from vol vars
+            unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+            velocityGradients_[elementID] = DimMatrix(0.1);
+            // TODO call kinematic viscosity value from vol vars
+            kinematicViscosity_[elementID] = 15e-6;
+        }
     }
 
     /*!
@@ -188,6 +205,8 @@ public:
 public:
     mutable std::vector<unsigned int> wallElementIDs_;
     mutable std::vector<Scalar> wallDistances_;
+    mutable std::vector<DimMatrix> velocityGradients_;
+    mutable std::vector<Scalar> kinematicViscosity_;
 
 private:
     //! Returns the implementation of the problem (i.e. static polymorphism)

dumux/freeflow/rans/volumevariables.hh

@@ -64,7 +64,8 @@ class RANSVolumeVariablesImplementation<TypeTag, false>
     using Element = typename GridView::template Codim<0>::Entity;
 
     enum { dimWorld = GridView::dimensionworld };
-    using FieldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+    using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+    using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 
     static const int defaultPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 
@@ -86,12 +87,20 @@ public:
                 const Element &element,
                 const SubControlVolume& scv)
     {
+        using std::sqrt;
         ParentType::update(elemSol, problem, element, scv);
 
+        asImp_().calculateVelocity(elemSol, problem, element, scv);
+        asImp_().calculateVelocityGradients(elemSol, problem, element, scv);
+
+        // calculate characteristic properties of the turbulent flow
         elementID_ = problem.fvGridGeometry().elementMapper().index(element);
+        wallElementID_ = problem.wallElementIDs_[elementID_];
         wallDistance_ = problem.wallDistances_[elementID_];
-
-        asImp_().calculateVelocityGradients(elemSol, problem, element, scv);
+        Scalar uStar = sqrt(problem.kinematicViscosity_[wallElementID_]
+                       * problem.velocityGradients_[wallElementID_][0][1]); // TODO: flow and wallnormalaxis
+        yPlus_ = wallDistance_ * uStar / asImp_().kinematicViscosity();
+        uPlus_ = velocity_[0] / uStar; // TODO: flow and wallnormalaxis
 
         // calculate the eddy viscosity based on the implemented RANS model
         asImp_().calculateEddyViscosity(elemSol, problem, element, scv);
@@ -99,7 +108,23 @@ public:
 
 
     /*!
-     * \brief Calculate the velocity gradients at the cell center
+     * \brief Calculate the velocity vector at the cell center
+     *
+     * \param elemSol A vector containing all primary variables connected to the element
+     * \param problem The object specifying the problem which ought to
+     *                be simulated
+     * \param element An element which contains part of the control volume
+     * \param scv The sub-control volume
+     */
+    void calculateVelocity(const ElementSolutionVector &elemSol,
+                           const Problem &problem,
+                           const Element &element,
+                           const SubControlVolume& scv)
+    { velocity_ = DimVector(1.0); /* TODO */ }
+
+
+    /*!
+     * \brief Calculate the velocity gradient tensor at the cell center
      *
      * \param elemSol A vector containing all primary variables connected to the element
      * \param problem The object specifying the problem which ought to
@@ -111,7 +136,7 @@ public:
                                     const Problem &problem,
                                     const Element &element,
                                     const SubControlVolume& scv)
-    { velocityGradients_ = FieldMatrix(1.0); }
+    { velocityGradients_ = DimMatrix(1.0); /* TODO */ }
 
 
     /*!
@@ -136,7 +161,7 @@ public:
     /*!
      * \brief Return the velocity gradients \f$\mathrm{[1/s]}\f$ at the control volume center.
      */
-    FieldMatrix velocityGradients() const
+    DimMatrix velocityGradients() const
     { return velocityGradients_; }
 
     /*!
@@ -145,6 +170,18 @@ public:
     Scalar wallDistance() const
     { return wallDistance_; }
 
+    /*!
+     * \brief Return the dimensionless wall distance \f$\mathrm{[-]}\f$.
+     */
+    Scalar yPlus() const
+    { return yPlus_; }
+
+    /*!
+     * \brief Return the dimensionless velocity \f$\mathrm{[-]}\f$.
+     */
+    Scalar uPlus() const
+    { return uPlus_; }
+
     /*!
      * \brief Set the values of the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$ within the
      *        control volume.
@@ -163,8 +200,15 @@ public:
      * \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
      *        control volume.
      */
-    Scalar effectiveViscosity(int phaseIdx = 0) const
-    { return asImp_().viscosity(defaultPhaseIdx) + dynamicEddyViscosity(); }
+    Scalar effectiveViscosity() const
+    { return asImp_().viscosity() + dynamicEddyViscosity(); }
+
+    /*!
+     * \brief Return the kinematic eddy viscosity \f$\mathrm{[m^2/s]}\f$ of the fluid within the
+     *        control volume.
+     */
+    Scalar kinematicViscosity() const
+    { return asImp_().viscosity() / asImp_().density(); }
 
 private:
     //! Returns the implementation of the problem (i.e. static polymorphism)
@@ -176,11 +220,14 @@ private:
     { return *static_cast<const Implementation *>(this); }
 
 protected:
-    FieldMatrix velocityGradients_;
+    DimVector velocity_;
+    DimMatrix velocityGradients_;
     Scalar dynamicEddyViscosity_;
     unsigned int elementID_;
     unsigned int wallElementID_;
     Scalar wallDistance_;
+    Scalar yPlus_;
+    Scalar uPlus_;
 };
 
 /*!
@@ -248,9 +295,9 @@ public:
      * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
      *        of the fluid-flow in the sub-control volume.
      */
-    Scalar effectiveThermalConductivity(int phaseIdx = 0) const
+    Scalar effectiveThermalConductivity() const
     {
-        return FluidSystem::thermalConductivity(this->fluidState_, defaultPhaseIdx)
+        return FluidSystem::thermalConductivity(this->fluidState_)
                + eddyThermalConductivity();
     }
 

dumux/freeflow/rans/vtkoutputfields.hh

@@ -64,6 +64,8 @@ public:
         vtk.addVolumeVariable([](const VolumeVariables& v){ return v.viscosity() / v.density(); }, "nu [m^2/s]");
         vtk.addVolumeVariable([](const VolumeVariables& v){ return v.dynamicEddyViscosity() / v.density(); }, "nu_t [m^2/s]");
         vtk.addVolumeVariable([](const VolumeVariables& v){ return v.wallDistance(); }, "l_w [m]");
+        vtk.addVolumeVariable([](const VolumeVariables& v){ return v.yPlus(); }, "y^+ [-]");
+        vtk.addVolumeVariable([](const VolumeVariables& v){ return v.uPlus(); }, "u^+ [-]");
 
         // add discretization-specific fields
         additionalOutput_(vtk, discMethodTag<GET_PROP_VALUE(TypeTag, DiscretizationMethod)>{});

dumux/freeflow/rans/zeroeq/CMakeLists.txt

-# #install headers
-# install(FILES
-# fluxoverplane.hh
-# fluxvariables.hh
-# localresidual.hh
-# DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/zeroeq)
+#install headers
+install(FILES
+indices.hh
+model.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/rans/zeroeq)

dumux/freeflow/rans/zeroeq/indices.hh

@@ -33,6 +33,8 @@ namespace Dumux
  * \ingroup ZeroEqModel
  * \brief The common indices for the isothermal ZeroEq model.
  *
+ * \tparam dimension The dimension of the problem
+ * \tparam numEquations the number of model equations
  * \tparam PVOffset The first index in a primary variable vector.
  */
 template <int dimension, int numEquations, int PVOffset = 0>
@@ -41,6 +43,7 @@ struct ZeroEqIndices
 {
     static constexpr int noEddyViscosityModel = 0;
     static constexpr int prandtl = 1;
+    static constexpr int modifiedVanDriest = 2;
 };
 
 // \}

dumux/freeflow/rans/zeroeq/volumevariables.hh

@@ -104,14 +104,32 @@ public:
                                 const Element &element,
                                 const SubControlVolume& scv)
     {
+        using std::exp;
         Scalar kinematicEddyViscosity = 0.0;
-        if (eddyViscosityModel_ == Indices::prandtl)
+        const Scalar karmanConstant = 0.41; // TODO make karman constant a property
+        if (eddyViscosityModel_ == Indices::noEddyViscosityModel)
         {
-            Scalar mixingLength = 0.41 * asImp_().wallDistance_;
-            Scalar velGrad = asImp_().velocityGradients()[0][1];
+            // kinematicEddyViscosity = 0.0
+        }
+        else if (eddyViscosityModel_ == Indices::prandtl)
+        {
+            Scalar mixingLength = karmanConstant * asImp_().wallDistance();
+            Scalar velGrad = asImp_().velocityGradients()[0][1]; // TODO: flow and wallnormalaxis
             kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
         }
-        asImp_().setDynamicEddyViscosity(kinematicEddyViscosity * asImp_().density(defaultPhaseIdx));
+        else if (eddyViscosityModel_ == Indices::modifiedVanDriest)
+        {
+            Scalar mixingLength = karmanConstant * asImp_().wallDistance()
+                                  * (1.0 - exp(-asImp_().yPlus() / 26.0));
+            Scalar velGrad = asImp_().velocityGradients()[0][1]; // TODO: flow and wallnormalaxis
+            kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
+        }
+        else
+        {
+            DUNE_THROW(Dune::NotImplemented,
+                       "This eddy viscosity model is not implemented: " << eddyViscosityModel_);
+        }
+        asImp_().setDynamicEddyViscosity(kinematicEddyViscosity * asImp_().density());
     }
 
 private:

test/freeflow/rans/test_pipe_laufer.cc

@@ -177,34 +177,6 @@ int main(int argc, char** argv) try
     using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
     NewtonSolver nonLinearSolver(assembler, linearSolver);
 
-    // set up two planes over which fluxes are calculated
-    FluxOverPlane<TypeTag> flux(*assembler, x);
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-    using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
-
-    const Scalar xMin = fvGridGeometry->bBoxMin()[0];
-    const Scalar xMax = fvGridGeometry->bBoxMax()[0];
-    const Scalar yMin = fvGridGeometry->bBoxMin()[1];
-    const Scalar yMax = fvGridGeometry->bBoxMax()[1];
-
-    // The first plane shall be placed at the middle of the channel.
-    // If we have an odd number of cells in x-direction, there would not be any cell faces
-    // at the postion of the plane (which is required for the flux calculation).
-    // In this case, we add half a cell-width to the x-position in order to make sure that
-    // the cell faces lie on the plane. This assumes a regular cartesian grid.
-    const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin);
-    const int numCellsX = getParam<std::vector<int>>("Grid.Cells")[0];
-    const Scalar offsetX = (numCellsX % 2 == 0) ? 0.0 : 0.5*((xMax - xMin) / numCellsX);
-
-    const auto p0middle = GlobalPosition{planePosMiddleX + offsetX, yMin};
-    const auto p1middle = GlobalPosition{planePosMiddleX + offsetX, yMax};
-    flux.addPlane("middle", p0middle, p1middle);
-
-    // The second plane is placed at the outlet of the channel.
-    const auto p0outlet = GlobalPosition{xMax, yMin};
-    const auto p1outlet = GlobalPosition{xMax, yMax};
-    flux.addPlane("outlet", p0outlet, p1outlet);
-
     // time loop
     timeLoop->start(); do
     {
@@ -218,6 +190,9 @@ int main(int argc, char** argv) try
         xOld = x;
         gridVariables->advanceTimeStep();
 
+        // update wall properties
+        problem->updateDynamicWallProperties();
+
         // advance to the time loop to the next step
         timeLoop->advanceTimeStep();