diff --git a/dumux/freeflow/rans/CMakeLists.txt b/dumux/freeflow/rans/CMakeLists.txt
index 74c1bf6c2ba729fcae9c2fe618fd9df57b69276c..b1ccda770bea551460d825be7284ae307a75d379 100644
--- a/dumux/freeflow/rans/CMakeLists.txt
+++ b/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
diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index bfc820c7ec34fb174a116d7e2ec51ff1aa520e2a..e645a7d930c49b288655fd0f957531900c3b72ee 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/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)
diff --git a/dumux/freeflow/rans/volumevariables.hh b/dumux/freeflow/rans/volumevariables.hh
index b0f668dbb871f0f0a47c01729bcca4230ebe5164..6decec4510a01a9b5db259dec3ab9add51a5dc65 100644
--- a/dumux/freeflow/rans/volumevariables.hh
+++ b/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();
}
diff --git a/dumux/freeflow/rans/vtkoutputfields.hh b/dumux/freeflow/rans/vtkoutputfields.hh
index 37de517e6129fcc9a12ea23d66b5b11c5eef586e..be6c9803dca16e26fbf72a28f932cf38bcc8112e 100644
--- a/dumux/freeflow/rans/vtkoutputfields.hh
+++ b/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)>{});
diff --git a/dumux/freeflow/rans/zeroeq/CMakeLists.txt b/dumux/freeflow/rans/zeroeq/CMakeLists.txt
index eb0484b81428c67c85c01379f04ca90add176cc4..811156a459b9ad0ac3fa9a641e2162c6cfa3aa75 100644
--- a/dumux/freeflow/rans/zeroeq/CMakeLists.txt
+++ b/dumux/freeflow/rans/zeroeq/CMakeLists.txt
@@ -1,6 +1,6 @@
-# #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)
diff --git a/dumux/freeflow/rans/zeroeq/indices.hh b/dumux/freeflow/rans/zeroeq/indices.hh
index 8b6c2e0145922bbf18dee10665922d7de4a8c86d..e37c261077d14dfafbfe46d2feffc6c02c6b897d 100644
--- a/dumux/freeflow/rans/zeroeq/indices.hh
+++ b/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;
};
// \}
diff --git a/dumux/freeflow/rans/zeroeq/volumevariables.hh b/dumux/freeflow/rans/zeroeq/volumevariables.hh
index 4de733d7d10c63f7d7375571091805c748a38d32..c5121828893dd3a8082fcef02edc65e45edbacdf 100644
--- a/dumux/freeflow/rans/zeroeq/volumevariables.hh
+++ b/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:
diff --git a/test/freeflow/rans/test_pipe_laufer.cc b/test/freeflow/rans/test_pipe_laufer.cc
index 25e544171ad3af93a5ca48d75bdfb9e416569f65..d1a5f7867e3c376de58db7a274b8e4a8001cb346 100644
--- a/test/freeflow/rans/test_pipe_laufer.cc
+++ b/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();