From 6c99d7f1cbb097076e255a9ea44854409b526bbb Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Mon, 19 Mar 2018 10:39:49 +0100
Subject: [PATCH] [rans][zeroeq] Implement roughness and various changes
* do NOT use updateStaticWallProperties() in problem constructor
* correct calculation of maximum velocity
* add Laufer velocity profile
* update wall variable names
---
dumux/freeflow/rans/problem.hh | 73 ++++++++------
dumux/freeflow/rans/volumevariables.hh | 35 +++++--
dumux/freeflow/rans/vtkoutputfields.hh | 1 +
dumux/freeflow/rans/zeroeq/problem.hh | 95 +++++++++++++++----
dumux/freeflow/rans/zeroeq/volumevariables.hh | 28 +++++-
test/freeflow/rans/CMakeLists.txt | 2 +-
test/freeflow/rans/laufer_re50000.csv | 42 ++++++++
test/freeflow/rans/laufer_re50000_u+y+.csv | 2 +-
test/freeflow/rans/pipelauferproblem.hh | 7 ++
test/freeflow/rans/test_pipe_laufer.cc | 29 ++++--
test/freeflow/rans/test_pipe_laufer.input | 4 +-
test/references/pipe_laufer_zeroeq.vtu | 26 ++++-
12 files changed, 270 insertions(+), 74 deletions(-)
create mode 100644 test/freeflow/rans/laufer_re50000.csv
diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index 3eb655df03..ef1b32a83e 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/dumux/freeflow/rans/problem.hh
@@ -86,9 +86,7 @@ public:
//! The constructor sets the gravity, if desired by the user.
RANSProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
- {
- updateStaticWallProperties();
- }
+ { }
/*!
* \brief Update the static (solution independent) relations to the walls
@@ -102,17 +100,16 @@ public:
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(), std::numeric_limits<Scalar>::max());
- neighborIDs_.resize(this->fvGridGeometry().elementMapper().size());
- cellCenters_.resize(this->fvGridGeometry().elementMapper().size(), GlobalPosition(0.0));
+ wallElementID_.resize(this->fvGridGeometry().elementMapper().size());
+ wallDistance_.resize(this->fvGridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max());
+ neighborID_.resize(this->fvGridGeometry().elementMapper().size());
+ cellCenter_.resize(this->fvGridGeometry().elementMapper().size(), GlobalPosition(0.0));
velocity_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(0.0));
- velocityMaximum_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(0.0));
- velocityMinimum_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
velocityGradients_.resize(this->fvGridGeometry().elementMapper().size(), DimMatrix(0.0));
flowNormalAxis_.resize(this->fvGridGeometry().elementMapper().size(), 0);
wallNormalAxis_.resize(this->fvGridGeometry().elementMapper().size(), 1);
kinematicViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+ sandGrainRoughness_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
// retrieve all wall intersections and corresponding elements
std::vector<unsigned int> wallElements;
@@ -142,7 +139,7 @@ public:
for (const auto& element : elements(gridView))
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
- cellCenters_[elementID] = element.geometry().center();
+ cellCenter_[elementID] = element.geometry().center();
for (unsigned int i = 0; i < wallPositions.size(); ++i)
{
static const int problemWallNormalAxis
@@ -158,13 +155,14 @@ public:
GlobalPosition global = element.geometry().center();
global -= wallPositions[i];
// second and argument ensures to use only aligned elements
- if (abs(global[searchAxis]) < wallDistances_[elementID]
+ if (abs(global[searchAxis]) < wallDistance_[elementID]
&& abs(global[searchAxis]) < global.two_norm() + 1e-8
&& abs(global[searchAxis]) > global.two_norm() - 1e-8)
{
- wallDistances_[elementID] = abs(global[searchAxis]);
- wallElementIDs_[elementID] = wallElements[i];
+ wallDistance_[elementID] = abs(global[searchAxis]);
+ wallElementID_[elementID] = wallElements[i];
wallNormalAxis_[elementID] = searchAxis;
+ sandGrainRoughness_[elementID] = asImp_().sandGrainRoughnessAtPos(wallPositions[i]);
}
}
}
@@ -175,8 +173,8 @@ public:
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
- neighborIDs_[elementID][dimIdx][0] = elementID;
- neighborIDs_[elementID][dimIdx][1] = elementID;
+ neighborID_[elementID][dimIdx][0] = elementID;
+ neighborID_[elementID][dimIdx][1] = elementID;
}
for (const auto& intersection : intersections(gridView, element))
@@ -187,15 +185,15 @@ public:
unsigned int neighborID = this->fvGridGeometry().elementMapper().index(intersection.outside());
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
- if (abs(cellCenters_[elementID][dimIdx] - cellCenters_[neighborID][dimIdx]) > 1e-8)
+ if (abs(cellCenter_[elementID][dimIdx] - cellCenter_[neighborID][dimIdx]) > 1e-8)
{
- if (cellCenters_[elementID][dimIdx] > cellCenters_[neighborID][dimIdx])
+ if (cellCenter_[elementID][dimIdx] > cellCenter_[neighborID][dimIdx])
{
- neighborIDs_[elementID][dimIdx][0] = neighborID;
+ neighborID_[elementID][dimIdx][0] = neighborID;
}
- if (cellCenters_[elementID][dimIdx] < cellCenters_[neighborID][dimIdx])
+ if (cellCenter_[elementID][dimIdx] < cellCenter_[neighborID][dimIdx])
{
- neighborIDs_[elementID][dimIdx][1] = neighborID;
+ neighborID_[elementID][dimIdx][1] = neighborID;
}
}
}
@@ -222,6 +220,10 @@ public:
static const int flowNormalAxis
= getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.FlowNormalAxis", -1);
+ // re-initialize min and max values
+ velocityMaximum_.assign(this->fvGridGeometry().elementMapper().size(), DimVector(0.0));
+ velocityMinimum_.assign(this->fvGridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
+
// calculate cell-center-averaged velocities
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
@@ -245,7 +247,7 @@ public:
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
- unsigned int wallElementID = wallElementIDs_[elementID];
+ unsigned int wallElementID = wallElementID_[elementID];
Scalar maxVelocity = 0.0;
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
@@ -253,10 +255,10 @@ public:
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
{
velocityGradients_[elementID][velIdx][dimIdx]
- = (velocity_[neighborIDs_[elementID][dimIdx][1]][velIdx]
- - velocity_[neighborIDs_[elementID][dimIdx][0]][velIdx])
- / (cellCenters_[neighborIDs_[elementID][dimIdx][1]][dimIdx]
- - cellCenters_[neighborIDs_[elementID][dimIdx][0]][dimIdx]);
+ = (velocity_[neighborID_[elementID][dimIdx][1]][velIdx]
+ - velocity_[neighborID_[elementID][dimIdx][0]][velIdx])
+ / (cellCenter_[neighborID_[elementID][dimIdx][1]][dimIdx]
+ - cellCenter_[neighborID_[elementID][dimIdx][0]][dimIdx]);
}
if (abs(velocity_[elementID][dimIdx]) > abs(velocityMaximum_[wallElementID][dimIdx]))
@@ -302,7 +304,7 @@ public:
/*!
* \brief Returns whether a given point is on a wall
*
- * \param globalPos The position in global coordinates where the temperature should be specified.
+ * \param globalPos The position in global coordinates.
*/
bool isOnWall(const GlobalPosition &globalPos) const
{
@@ -311,11 +313,21 @@ public:
"The problem does not provide an isOnWall() method.");
}
+ /*!
+ * \brief Returns the sand-grain roughness \f$\mathrm{[m]}\f$ at a given position
+ *
+ * \param globalPos The position in global coordinates.
+ */
+ Scalar sandGrainRoughnessAtPos(const GlobalPosition &globalPos) const
+ {
+ return 0.0;
+ }
+
public:
- std::vector<unsigned int> wallElementIDs_;
- std::vector<Scalar> wallDistances_;
- std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIDs_;
- std::vector<GlobalPosition> cellCenters_;
+ std::vector<unsigned int> wallElementID_;
+ std::vector<Scalar> wallDistance_;
+ std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborID_;
+ std::vector<GlobalPosition> cellCenter_;
std::vector<DimVector> velocity_;
std::vector<DimVector> velocityMaximum_;
std::vector<DimVector> velocityMinimum_;
@@ -323,6 +335,7 @@ public:
std::vector<unsigned int> wallNormalAxis_;
std::vector<unsigned int> flowNormalAxis_;
std::vector<Scalar> kinematicViscosity_;
+ std::vector<Scalar> sandGrainRoughness_;
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 828b79a5a6..1bc2bfdf59 100644
--- a/dumux/freeflow/rans/volumevariables.hh
+++ b/dumux/freeflow/rans/volumevariables.hh
@@ -101,25 +101,38 @@ public:
// calculate characteristic properties of the turbulent flow
elementID_ = problem.fvGridGeometry().elementMapper().index(element);
- wallElementID_ = problem.wallElementIDs_[elementID_];
- wallDistance_ = problem.wallDistances_[elementID_];
+ wallElementID_ = problem.wallElementID_[elementID_];
+ wallDistance_ = problem.wallDistance_[elementID_];
velocity_ = problem.velocity_[elementID_];
+ velocityMaximum_ = problem.velocityMaximum_[wallElementID_];
velocityGradients_ = problem.velocityGradients_[elementID_];
unsigned int flowNormalAxis = problem.flowNormalAxis_[elementID_];
unsigned int wallNormalAxis = problem.wallNormalAxis_[elementID_];
- Scalar uStar = sqrt(problem.kinematicViscosity_[wallElementID_]
- * abs(problem.velocityGradients_[wallElementID_][flowNormalAxis][wallNormalAxis]));
- yPlus_ = wallDistance_ * uStar / asImp_().kinematicViscosity();
- yPlus_ = max(yPlus_, 1e-10); // zero values lead to numerical problems in some turbulence models
- uPlus_ = velocity_[flowNormalAxis] / max(uStar, 1e-10);
+ uStar_ = sqrt(problem.kinematicViscosity_[wallElementID_]
+ * abs(problem.velocityGradients_[wallElementID_][flowNormalAxis][wallNormalAxis]));
+ uStar_ = max(uStar_, 1e-10); // zero values lead to numerical problems in some turbulence models
+ yPlus_ = wallDistance_ * uStar_ / asImp_().kinematicViscosity();
+ uPlus_ = velocity_[flowNormalAxis] / uStar_;
};
+ /*!
+ * \brief Return the element ID of the control volume.
+ */
+ unsigned int elementID() const
+ { return elementID_; }
+
/*!
* \brief Return the velocity vector \f$\mathrm{[m/s]}\f$ at the control volume center.
*/
DimVector velocity() const
{ return velocity_; }
+ /*!
+ * \brief Return the maximum velocity vector \f$\mathrm{[m/s]}\f$ of the wall segment.
+ */
+ DimVector velocityMaximum() const
+ { return velocityMaximum_; }
+
/*!
* \brief Return the velocity gradients \f$\mathrm{[1/s]}\f$ at the control volume center.
*/
@@ -132,6 +145,12 @@ public:
Scalar wallDistance() const
{ return wallDistance_; }
+ /*!
+ * \brief Return the wall friction velocity \f$\mathrm{[m/s]}\f$
+ */
+ Scalar uStar() const
+ { return uStar_; }
+
/*!
* \brief Return the dimensionless wall distance \f$\mathrm{[-]}\f$.
*/
@@ -183,11 +202,13 @@ private:
protected:
DimVector velocity_;
+ DimVector velocityMaximum_;
DimMatrix velocityGradients_;
Scalar dynamicEddyViscosity_;
unsigned int elementID_;
unsigned int wallElementID_;
Scalar wallDistance_;
+ Scalar uStar_;
Scalar yPlus_;
Scalar uPlus_;
};
diff --git a/dumux/freeflow/rans/vtkoutputfields.hh b/dumux/freeflow/rans/vtkoutputfields.hh
index ff84b83d83..b1a1d3d0a6 100644
--- a/dumux/freeflow/rans/vtkoutputfields.hh
+++ b/dumux/freeflow/rans/vtkoutputfields.hh
@@ -45,6 +45,7 @@ public:
{
NavierStokesVtkOutputFields<FVGridGeometry>::init(vtk);
+ vtk.addVolumeVariable([](const auto& v){ return v.velocity()[0] / v.velocityMaximum()[0]; }, "v_x/v_x,max");
vtk.addVolumeVariable([](const auto& v){ return v.velocityGradients()[0]; }, "dv_x/dx_");
if (dim > 1)
vtk.addVolumeVariable([](const auto& v){ return v.velocityGradients()[1]; }, "dv_y/dx_");
diff --git a/dumux/freeflow/rans/zeroeq/problem.hh b/dumux/freeflow/rans/zeroeq/problem.hh
index 2c0841d90c..53dccbe458 100644
--- a/dumux/freeflow/rans/zeroeq/problem.hh
+++ b/dumux/freeflow/rans/zeroeq/problem.hh
@@ -82,23 +82,24 @@ public:
//! The constructor sets the gravity, if desired by the user.
ZeroEqProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
- {
- updateStaticWallProperties();
- }
+ { }
/*!
* \brief Correct size of the static (solution independent) wall variables
*/
void updateStaticWallProperties()
{
+ ParentType::updateStaticWallProperties();
+
// update size and initial values of the global vectors
kinematicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
+ additionalRoughnessLength_.resize(this->fvGridGeometry().elementMapper().size(), 0.0);
}
/*!
* \brief Update the dynamic (solution dependent) relations to the walls
*
- * This calculate all necessary information for the Baldwin-Lomax turbulence model
+ * This calculates the roughness related properties
*
* \param curSol The solution vector.
*/
@@ -106,16 +107,67 @@ public:
{
ParentType::updateDynamicWallProperties(curSol);
+ static const int eddyViscosityModel
+ = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.EddyViscosityModel");
+
+ // calculate additional roughness
+ bool printedRangeWarning = false;
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ using std::sqrt;
+ using std::exp;
+
+ const int dofIdx = scv.dofIndex();
+ CellCenterPrimaryVariables priVars(curSol[cellCenterIdx][dofIdx]);
+ auto elemSol = elementSolution<SolutionVector, FVGridGeometry>(std::move(priVars));
+ VolumeVariables volVars;
+ volVars.update(elemSol, asImp_(), element, scv);
+
+ Scalar ksPlus = this->sandGrainRoughness_[elementID] * volVars.uStar() / volVars.kinematicViscosity();
+ if (ksPlus > 0 && eddyViscosityModel == Indices::baldwinLomax)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Roughness is not implemented for the Baldwin-Lomax model.");
+ }
+ if (ksPlus > 2000.)
+ {
+ std::cout << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << this->cellCenter_[this->wallElementID_[elementID]]
+ << " is not in the valid range (ksPlus < 2000),"
+ << " for high ksPlus values the roughness function reaches a turning point."<< std::endl;
+ DUNE_THROW(Dune::InvalidStateException, "Unphysical roughness behavior.");
+ }
+ else if (ksPlus > 0.0 && ksPlus < 4.535 && !printedRangeWarning)
+ {
+ Dune::dinfo << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << this->cellCenter_[this->wallElementID_[elementID]]
+ << " is not in the valid range (ksPlus > 4.535) and now set to 0.0"<< std::endl;
+ ksPlus = 0.0;
+ printedRangeWarning = true;
+ }
+ additionalRoughnessLength_[elementID] = 0.9 / (volVars.uStar() / volVars.kinematicViscosity())
+ * (sqrt(ksPlus) - ksPlus * exp(-ksPlus / 6.0));
+ }
+ }
+
+ // update routine for specfic models
+ if (eddyViscosityModel == Indices::baldwinLomax)
+ updateBaldwinLomaxProperties();
+ }
+
+ /*!
+ * \brief Update the relations and coefficients for the Baldwin-Lomax turbulence model
+ */
+ void updateBaldwinLomaxProperties()
+ {
std::vector<Scalar> kinematicEddyViscosityInner(this->fvGridGeometry().elementMapper().size(), 0.0);
std::vector<Scalar> kinematicEddyViscosityOuter(this->fvGridGeometry().elementMapper().size(), 0.0);
std::vector<Scalar> kinematicEddyViscosityDifference(this->fvGridGeometry().elementMapper().size(), 0.0);
std::vector<Scalar> switchingPosition(this->fvGridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max());
- static const int eddyViscosityModel
- = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "RANS.EddyViscosityModel");
- if (eddyViscosityModel != Indices::baldwinLomax)
- return;
-
using std::abs;
using std::exp;
using std::min;
@@ -139,8 +191,8 @@ public:
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
- unsigned int wallElementID = this->wallElementIDs_[elementID];
- Scalar wallDistance = this->wallDistances_[elementID];
+ unsigned int wallElementID = this->wallElementID_[elementID];
+ Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID];
unsigned int flowNormalAxis = this->flowNormalAxis_[elementID];
unsigned int wallNormalAxis = this->wallNormalAxis_[elementID];
@@ -164,8 +216,8 @@ public:
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
- unsigned int wallElementID = this->wallElementIDs_[elementID];
- Scalar wallDistance = this->wallDistances_[elementID];
+ unsigned int wallElementID = this->wallElementID_[elementID];
+ Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID];
Scalar maxVelocityNorm = 0.0;
Scalar minVelocityNorm = 0.0;
@@ -191,8 +243,8 @@ public:
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
- unsigned int wallElementID = this->wallElementIDs_[elementID];
- Scalar wallDistance = this->wallDistances_[elementID];
+ unsigned int wallElementID = this->wallElementID_[elementID];
+ Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID];
// checks if sign switches, by multiplication
Scalar check = kinematicEddyViscosityDifference[wallElementID] * kinematicEddyViscosityDifference[elementID];
@@ -207,17 +259,20 @@ public:
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
- unsigned int wallElementID = this->wallElementIDs_[elementID];
- Scalar wallDistance = this->wallDistances_[elementID];
+ unsigned int wallElementID = this->wallElementID_[elementID];
+ Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID];
- kinematicEddyViscosity_[elementID] = (wallDistance >= switchingPosition[wallElementID])
- ? kinematicEddyViscosityOuter[elementID]
- : kinematicEddyViscosityInner[elementID];
+ kinematicEddyViscosity_[elementID] = kinematicEddyViscosityInner[elementID];
+ if (wallDistance >= switchingPosition[wallElementID])
+ {
+ kinematicEddyViscosity_[elementID] = kinematicEddyViscosityOuter[elementID];
+ }
}
}
public:
std::vector<Scalar> kinematicEddyViscosity_;
+ std::vector<Scalar> additionalRoughnessLength_;
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
diff --git a/dumux/freeflow/rans/zeroeq/volumevariables.hh b/dumux/freeflow/rans/zeroeq/volumevariables.hh
index f933af541d..0e80250f2c 100644
--- a/dumux/freeflow/rans/zeroeq/volumevariables.hh
+++ b/dumux/freeflow/rans/zeroeq/volumevariables.hh
@@ -83,10 +83,12 @@ public:
const SubControlVolume& scv)
{
ParentType::update(elemSol, problem, element, scv);
+ additionalRoughnessLength_ = problem.additionalRoughnessLength_[this->elementID()];
+ yPlusRough_ = wallDistanceRough() * this->uStar() / this->kinematicViscosity();
+
calculateEddyViscosity(elemSol, problem, element, scv);
}
-
/*!
* \brief Calculate and set the dynamic eddy viscosity.
*
@@ -121,14 +123,14 @@ public:
}
else if (eddyViscosityModel == Indices::prandtl)
{
- Scalar mixingLength = karmanConstant * asImp_().wallDistance();
+ Scalar mixingLength = karmanConstant * asImp_().wallDistanceRough();
kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
}
else if (eddyViscosityModel == Indices::modifiedVanDriest)
{
- Scalar mixingLength = karmanConstant * asImp_().wallDistance()
- * (1.0 - exp(-asImp_().yPlus() / 26.0))
- / sqrt(1.0 - exp(-0.26 * asImp_().yPlus()));
+ Scalar mixingLength = karmanConstant * asImp_().wallDistanceRough()
+ * (1.0 - exp(-asImp_().yPlusRough() / 26.0))
+ / sqrt(1.0 - exp(-0.26 * asImp_().yPlusRough()));
kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
}
else if (eddyViscosityModel == Indices::baldwinLomax)
@@ -143,6 +145,18 @@ public:
asImp_().setDynamicEddyViscosity(kinematicEddyViscosity * asImp_().density());
}
+ /*!
+ * \brief Return the wall distance \f$\mathrm{[m]}\f$ including an additional roughness length
+ */
+ Scalar wallDistanceRough() const
+ { return this->wallDistance() + additionalRoughnessLength_; }
+
+ /*!
+ * \brief Return the dimensionless wall distance \f$\mathrm{[-]}\f$ including an additional roughness length
+ */
+ Scalar yPlusRough() const
+ { return yPlusRough_; }
+
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
@@ -151,6 +165,10 @@ private:
//! \copydoc asImp_()
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
+
+protected:
+ Scalar additionalRoughnessLength_;
+ Scalar yPlusRough_;
};
/*!
diff --git a/test/freeflow/rans/CMakeLists.txt b/test/freeflow/rans/CMakeLists.txt
index 0b3b205ce7..a9e36ffdb4 100644
--- a/test/freeflow/rans/CMakeLists.txt
+++ b/test/freeflow/rans/CMakeLists.txt
@@ -1,5 +1,5 @@
add_input_file_links()
-dune_symlink_to_source_files(FILES laufer_re50000_u+y+.csv)
+dune_symlink_to_source_files(FILES laufer_re50000.csv laufer_re50000_u+y+.csv)
dune_add_test(NAME test_pipe_laufer
SOURCES test_pipe_laufer.cc
diff --git a/test/freeflow/rans/laufer_re50000.csv b/test/freeflow/rans/laufer_re50000.csv
new file mode 100644
index 0000000000..d695d674e3
--- /dev/null
+++ b/test/freeflow/rans/laufer_re50000.csv
@@ -0,0 +1,42 @@
+###########
+# Velocity profile for Re=500000
+###########
+# Data taken from:
+# Laufer, J.
+# The structure of turbulence in fully developed pipe flow
+# NACA Report, 1954, 1174, 417-434
+###########
+#y/d,u/u_max
+0.000,0.000
+0.006,0.333
+0.008,0.480
+0.014,0.549
+0.019,0.617
+0.033,0.648
+0.048,0.699
+0.078,0.740
+0.107,0.792
+0.141,0.831
+0.195,0.866
+0.244,0.902
+0.295,0.932
+0.348,0.959
+0.401,0.975
+0.452,0.988
+0.500,1.000
+0.549,0.988
+0.599,0.975
+0.652,0.959
+0.706,0.932
+0.757,0.902
+0.806,0.866
+0.860,0.831
+0.893,0.792
+0.923,0.740
+0.953,0.699
+0.968,0.648
+0.981,0.617
+0.987,0.549
+0.992,0.480
+0.995,0.333
+1.000,0.000
diff --git a/test/freeflow/rans/laufer_re50000_u+y+.csv b/test/freeflow/rans/laufer_re50000_u+y+.csv
index 547f79ab0c..53788c4661 100644
--- a/test/freeflow/rans/laufer_re50000_u+y+.csv
+++ b/test/freeflow/rans/laufer_re50000_u+y+.csv
@@ -1,5 +1,5 @@
###########
-# Dimensionless velocity profile for Re=50,000
+# Dimensionless velocity profile for Re=50000
###########
# Original source:
# Laufer, J.
diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh
index a9b0e0c68b..ea8e99b495 100644
--- a/test/freeflow/rans/pipelauferproblem.hh
+++ b/test/freeflow/rans/pipelauferproblem.hh
@@ -105,6 +105,7 @@ public:
: ParentType(fvGridGeometry), eps_(1e-6)
{
inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
+ sandGrainRoughness_ = getParam<Scalar>("Problem.SandGrainRoughness", 0.0);
}
/*!
@@ -119,6 +120,11 @@ public:
|| globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - localEps_;
}
+ Scalar sandGrainRoughnessAtPos(const GlobalPosition &globalPos) const
+ {
+ return sandGrainRoughness_;
+ }
+
bool shouldWriteRestartFile() const
{
return false;
@@ -228,6 +234,7 @@ private:
Scalar eps_;
Scalar inletVelocity_;
+ Scalar sandGrainRoughness_;
TimeLoopPtr timeLoop_;
};
} //end namespace
diff --git a/test/freeflow/rans/test_pipe_laufer.cc b/test/freeflow/rans/test_pipe_laufer.cc
index 261ff49ad4..a348a46983 100644
--- a/test/freeflow/rans/test_pipe_laufer.cc
+++ b/test/freeflow/rans/test_pipe_laufer.cc
@@ -142,6 +142,7 @@ int main(int argc, char** argv) try
x[cellCenterIdx].resize(numDofsCellCenter);
x[faceIdx].resize(numDofsFace);
problem->applyInitialSolution(x);
+ problem->updateStaticWallProperties();
problem->updateDynamicWallProperties(x);
auto xOld = x;
@@ -205,8 +206,9 @@ int main(int argc, char** argv) try
////////////////////////////////////////////////////////////
#if HAVE_PVPYTHON
- bool shouldPlot = getParam<bool>("Output.PlotVelocityProfiles", false);
- if (shouldPlot)
+ bool plotLawOfTheWall = getParam<bool>("Output.PlotLawOfTheWall", false);
+ bool plotVelocityProfile = getParam<bool>("Output.PlotVelocityProfile", false);
+ if (plotLawOfTheWall || plotVelocityProfile)
{
char fileName[255];
std::string fileNameFormat = "%s-%05d";
@@ -219,26 +221,37 @@ int main(int argc, char** argv) try
// execute the pvpython script
std::string command = std::string(PVPYTHON_EXECUTABLE) + " " + script
+ " -f " + vtuFileName
- + " -v 1"
+ + " -v 2"
+ " -r 10000";
syscom = command + " -p1 8.0 0.0 0.0"
+ " -p2 8.0 0.2469 0.0"
+ " -of " + std::string(fileName) + "\n";
+// std::cout << syscom << std::endl;
system(syscom.c_str());
-
char gnuplotFileName[255];
Dumux::GnuplotInterface<Scalar> gnuplot;
- sprintf(gnuplotFileName, fileNameFormat.c_str(), "velProfiles", timeLoop->timeStepIndex());
- gnuplot.setOpenPlotWindow(shouldPlot);
+ sprintf(gnuplotFileName, fileNameFormat.c_str(), "lawOfTheWall", timeLoop->timeStepIndex());
+ gnuplot.setOpenPlotWindow(plotLawOfTheWall);
gnuplot.setDatafileSeparator(',');
gnuplot.resetPlot();
gnuplot.setXlabel("y^+ [-]");
gnuplot.setYlabel("u_+ [-]");
gnuplot.setOption("set log x");
gnuplot.setOption("set xrange [1:3000]");
- gnuplot.addFileToPlot("laufer_re50000_u+y+.csv", "u 1:2 w p t 'Laufer1954a, Re=50000'");
- gnuplot.addFileToPlot(std::string(fileName) + ".csv", "u 10:11 w l");
+ gnuplot.addFileToPlot("laufer_re50000_u+y+.csv", "u 1:2 w p t 'Laufer 1954, Re=50000'");
+ gnuplot.addFileToPlot(std::string(fileName) + ".csv", "u 11:12 w l");
+ gnuplot.plot(std::string(gnuplotFileName));
+
+ sprintf(gnuplotFileName, fileNameFormat.c_str(), "velProfile", timeLoop->timeStepIndex());
+ gnuplot.resetAll();
+ gnuplot.setOpenPlotWindow(plotVelocityProfile);
+ gnuplot.setDatafileSeparator(',');
+ gnuplot.setXlabel("v_x/v_{x,max} [-]");
+ gnuplot.setYRange(0.0, 1.0);
+ gnuplot.setYlabel("y [-]");
+ gnuplot.addFileToPlot("laufer_re50000.csv", "u 2:1 w p t 'Laufer 1954, Re=50000'");
+ gnuplot.addFileToPlot(std::string(fileName) + ".csv", "u 5:($23/0.2456) w l");
gnuplot.plot(std::string(gnuplotFileName));
}
#endif
diff --git a/test/freeflow/rans/test_pipe_laufer.input b/test/freeflow/rans/test_pipe_laufer.input
index f507011c75..9acb019e21 100644
--- a/test/freeflow/rans/test_pipe_laufer.input
+++ b/test/freeflow/rans/test_pipe_laufer.input
@@ -11,12 +11,14 @@ Cells1 = 25 25
Grading1 = 1.2 -1.2
[Output]
-PlotVelocityProfiles = true
+PlotLawOfTheWall = true
+PlotVelocityProfile = true
[Problem]
Name = pipe_laufer_zeroeq
InletVelocity = 2.5 # [m/s]
EnableGravity = false
+SandGrainRoughness = 0.0 # [m] # not implemented for EddyViscosityModel = 3
[RANS]
EddyViscosityModel = 3
diff --git a/test/references/pipe_laufer_zeroeq.vtu b/test/references/pipe_laufer_zeroeq.vtu
index 4f3ada09cb..d1a4169806 100644
--- a/test/references/pipe_laufer_zeroeq.vtu
+++ b/test/references/pipe_laufer_zeroeq.vtu
@@ -2,7 +2,7 @@
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="256" NumberOfPoints="289">
- <CellData Scalars="p" Vectors="velocity_Air (m/s)">
+ <CellData Scalars="p" Vectors="dv_x/dx_">
<DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
100000 100000 100000 100000 100001 100001 100001 100001 100001 100001 100001 100001
@@ -27,6 +27,30 @@
100001 100001 100001 100001 100001 100001 100001 100001 100001 100000 100000 100000
100000 100000 100000 100000
</DataArray>
+ <DataArray type="Float32" Name="v_x/v_x,max" NumberOfComponents="1" format="ascii">
+ 0.824026 0.577694 0.494709 0.478182 0.462302 0.44743 0.440993 0.437014 0.431486 0.425143 0.42038 0.419126
+ 0.419458 0.419144 0.418159 0.417228 0.946449 0.837676 0.750475 0.702608 0.67421 0.658205 0.654306 0.651134
+ 0.643928 0.635299 0.628737 0.626977 0.626963 0.62594 0.62425 0.622795 0.984664 0.939124 0.877889 0.824945
+ 0.787385 0.762412 0.747611 0.736606 0.725543 0.714694 0.708117 0.708114 0.709018 0.7078 0.705826 0.704164
+ 0.997136 0.984401 0.956349 0.918926 0.883439 0.853441 0.829141 0.810303 0.794956 0.781691 0.772737 0.769462
+ 0.767816 0.765339 0.762618 0.760546 0.999713 0.997748 0.990327 0.975018 0.955061 0.933334 0.911037 0.890703
+ 0.873326 0.858481 0.846094 0.836777 0.83034 0.825476 0.821414 0.818646 1 0.999888 0.998978 0.995913
+ 0.990388 0.98199 0.970179 0.956788 0.943684 0.931463 0.920136 0.909932 0.901567 0.894987 0.88971 0.886277
+ 0.999988 1 1 0.999791 0.9993 0.998099 0.995474 0.991556 0.98692 0.981885 0.97636 0.97039
+ 0.964742 0.959863 0.955718 0.952956 0.999959 0.999947 0.999996 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 0.999959 0.999947 0.999996 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 0.999988 1 1 0.999791 0.9993 0.998099 0.995474 0.991556 0.98692 0.981885 0.97636 0.97039
+ 0.964742 0.959863 0.955718 0.952956 1 0.999888 0.998978 0.995913 0.990388 0.98199 0.970179 0.956788
+ 0.943684 0.931463 0.920136 0.909932 0.901567 0.894987 0.88971 0.886277 0.999713 0.997748 0.990327 0.975018
+ 0.955061 0.933334 0.911037 0.890703 0.873326 0.858481 0.846094 0.836777 0.83034 0.825476 0.821414 0.818646
+ 0.997136 0.984401 0.956349 0.918926 0.883439 0.853441 0.829141 0.810303 0.794956 0.781691 0.772737 0.769462
+ 0.767816 0.765339 0.762618 0.760546 0.984664 0.939124 0.877889 0.824945 0.787385 0.762412 0.747611 0.736606
+ 0.725543 0.714694 0.708117 0.708114 0.709018 0.7078 0.705826 0.704164 0.946449 0.837676 0.750475 0.702608
+ 0.67421 0.658205 0.654306 0.651134 0.643928 0.635299 0.628737 0.626977 0.626963 0.62594 0.62425 0.622795
+ 0.824026 0.577694 0.494709 0.478182 0.462302 0.44743 0.440993 0.437014 0.431486 0.425143 0.42038 0.419126
+ 0.419458 0.419144 0.418159 0.417228
+ </DataArray>
<DataArray type="Float32" Name="p_rel" NumberOfComponents="1" format="ascii">
1.30543 1.17114 1.0534 0.957444 0.867397 0.780484 0.700841 0.624126 0.545742 0.467086 0.389893 0.318209
0.249301 0.178674 0.107209 0.0356539 1.30533 1.17127 1.05352 0.957484 0.867414 0.780489 0.700836 0.624126
--
GitLab