From c0338acb1c6b76fe77ca23c106a7ebd9169c7214 Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Wed, 14 Mar 2018 14:43:42 +0100
Subject: [PATCH] [rans][zeroeq] First version of velocity gradients
---
dumux/freeflow/rans/problem.hh | 85 +++++++++++++++++--
dumux/freeflow/rans/volumevariables.hh | 7 +-
dumux/freeflow/rans/vtkoutputfields.hh | 2 +
dumux/freeflow/rans/zeroeq/volumevariables.hh | 8 +-
test/freeflow/rans/pipelauferproblem.hh | 16 +---
test/freeflow/rans/test_pipe_laufer.input | 3 +-
6 files changed, 98 insertions(+), 23 deletions(-)
diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index c188fd9f60..c49b403825 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/dumux/freeflow/rans/problem.hh
@@ -116,12 +116,15 @@ public:
// update size and initial values of the global vectors
wallElementIDs_.resize(this->fvGridGeometry().elementMapper().size());
wallDistances_.resize(this->fvGridGeometry().elementMapper().size());
+ neighborIDs_.resize(this->fvGridGeometry().elementMapper().size());
+ cellCenters_.resize(this->fvGridGeometry().elementMapper().size());
velocity_.resize(this->fvGridGeometry().elementMapper().size());
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();
+ cellCenters_[i] = GlobalPosition(0.0);
velocity_[i] = DimVector(0.0);
velocityGradients_[i] = DimMatrix(0.0);
kinematicViscosity_[i] = 0.0;
@@ -148,11 +151,12 @@ public:
// search for shortest distance to wall for each element
for (const auto& element : elements(gridView))
{
+ unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+ cellCenters_[elementID] = element.geometry().center();
for (unsigned int i = 0; i < wallPositions.size(); ++i)
{
GlobalPosition global = element.geometry().center();
global -= wallPositions[i];
- unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
if (global.two_norm() < wallDistances_[elementID])
{
wallDistances_[elementID] = global.two_norm();
@@ -160,6 +164,50 @@ public:
}
}
}
+
+ // search for neighbor IDs
+ for (const auto& element : elements(gridView))
+ {
+ unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+ std::array<std::array<Scalar, 2>, dim> distances;
+ for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+ {
+ neighborIDs_[elementID][dimIdx][0] = elementID;
+ neighborIDs_[elementID][dimIdx][1] = elementID;
+ distances[dimIdx][0] = std::numeric_limits<Scalar>::max();
+ distances[dimIdx][1] = -std::numeric_limits<Scalar>::max();
+ }
+ for (const auto& neighbor : elements(gridView))
+ {
+ unsigned int neighborID = this->fvGridGeometry().elementMapper().index(neighbor);
+ for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+ {
+ if (elementID == neighborID)
+ continue;
+// std::cout << elementID << " " << neighborID << std::endl;
+
+ for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+ {
+ Scalar distanceTemp = cellCenters_[elementID][dimIdx] - cellCenters_[neighborID][dimIdx];
+// std::cout << distanceTemp
+// << " " << distances[dimIdx][0]
+// << " " << distances[dimIdx][1]
+// << std::endl;
+ if (distanceTemp < distances[dimIdx][0] && distanceTemp > 1e-8)
+ {
+ neighborIDs_[elementID][dimIdx][0] = neighborID;
+ distances[dimIdx][0] = distanceTemp;
+ }
+
+ if (distanceTemp > distances[dimIdx][1] && distanceTemp < -1e-8)
+ {
+ neighborIDs_[elementID][dimIdx][1] = neighborID;
+ distances[dimIdx][1] = distanceTemp;
+ }
+ }
+ }
+ }
+ }
}
/*!
@@ -186,7 +234,29 @@ public:
}
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
velocity_[elementID][dimIdx] = velocityTemp[dimIdx] * 0.5; // faces are equidistant to cell center
+ }
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
+ for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+ {
+ 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]);
+// std::cout << " velocity_[1][velIdx] " << velocity_[neighborIDs_[elementID][dimIdx][1]][velIdx]
+// << " velocity_[0][velIdx] " << velocity_[neighborIDs_[elementID][dimIdx][0]][velIdx]
+// << " cellCenters_[1][velIdx] " << cellCenters_[neighborIDs_[elementID][dimIdx][1]][dimIdx]
+// << " cellCenters_[0][velIdx] " << cellCenters_[neighborIDs_[elementID][dimIdx][0]][dimIdx]
+// << " velocityGradients_[elementID][" << velIdx << "][" << dimIdx << "] " << velocityGradients_[elementID][velIdx][dimIdx];
+ }
+ }
+// std::cout << std::endl;
+ }
// // TODO: calculate velocity gradients
// for (auto&& scv : scvs(fvGeometry))
// {
@@ -218,12 +288,15 @@ public:
// velocityTemp[scvf.directionIndex()] += numericalSolutionFace;
// }
// }
- for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
- for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
- velocityGradients_[elementID][velIdx][dimIdx] = 0.12345;
+// for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
+// for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
+// velocityGradients_[elementID][velIdx][dimIdx] = 0.12345;
// // TODO calculate or call all secondary variables
// // TODO call kinematic viscosity value from vol vars
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ unsigned int elementID = this->fvGridGeometry().elementMapper().index(element);
kinematicViscosity_[elementID] = 15e-6;
}
}
@@ -231,7 +304,7 @@ public:
/*!
* \brief Returns whether a given point is on a wall
*/
- const bool isOnWall(GlobalPosition &globalPos) const
+ bool isOnWall(const GlobalPosition &globalPos) const
{
// Throw an exception if no walls are implemented
DUNE_THROW(Dune::InvalidStateException,
@@ -281,6 +354,8 @@ public:
public:
mutable std::vector<unsigned int> wallElementIDs_;
mutable std::vector<Scalar> wallDistances_;
+ mutable std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIDs_;
+ mutable std::vector<GlobalPosition> cellCenters_;
mutable std::vector<DimVector> velocity_;
mutable std::vector<DimMatrix> velocityGradients_;
mutable std::vector<Scalar> kinematicViscosity_;
diff --git a/dumux/freeflow/rans/volumevariables.hh b/dumux/freeflow/rans/volumevariables.hh
index 9851545ba0..d0df6196b9 100644
--- a/dumux/freeflow/rans/volumevariables.hh
+++ b/dumux/freeflow/rans/volumevariables.hh
@@ -87,6 +87,8 @@ public:
const Element &element,
const SubControlVolume& scv)
{
+ using std::abs;
+ using std::max;
using std::sqrt;
ParentType::update(elemSol, problem, element, scv);
@@ -95,10 +97,11 @@ public:
wallElementID_ = problem.wallElementIDs_[elementID_];
wallDistance_ = problem.wallDistances_[elementID_];
velocity_ = problem.velocity_[elementID_];
+ velocityGradients_ = problem.velocityGradients_[elementID_];
Scalar uStar = sqrt(problem.kinematicViscosity_[wallElementID_]
- * problem.velocityGradients_[wallElementID_][0][1]); // TODO: flow and wallnormalaxis
+ * abs(problem.velocityGradients_[wallElementID_][0][1])); // TODO: flow and wallnormalaxis
yPlus_ = wallDistance_ * uStar / asImp_().kinematicViscosity();
- uPlus_ = velocity_[0] / uStar; // TODO: flow and wallnormalaxis
+ uPlus_ = velocity_[0] / max(uStar, 1e-8); // TODO: flow and wallnormalaxis
// calculate the eddy viscosity based on the implemented RANS model
asImp_().calculateEddyViscosity(elemSol, problem, element, scv);
diff --git a/dumux/freeflow/rans/vtkoutputfields.hh b/dumux/freeflow/rans/vtkoutputfields.hh
index 40d5a8c472..631bb23006 100644
--- a/dumux/freeflow/rans/vtkoutputfields.hh
+++ b/dumux/freeflow/rans/vtkoutputfields.hh
@@ -52,6 +52,8 @@ public:
static void init(VtkOutputModule& vtk)
{
vtk.addVolumeVariable([](const auto& v){ return v.velocity()[0]; }, "v_x [m/s]");
+ vtk.addVolumeVariable([](const auto& v){ return v.velocityGradients()[0][0]; }, "dv_x/dx [m/s]");
+ vtk.addVolumeVariable([](const auto& v){ return v.velocityGradients()[0][1]; }, "dv_x/dy [m/s]");
vtk.addVolumeVariable([](const auto& v){ return v.pressure(); }, "p [Pa]");
vtk.addVolumeVariable([](const auto& v){ return v.pressure() - 1e5; }, "p_rel [Pa]");
vtk.addVolumeVariable([](const auto& v){ return v.density(); }, "rho [kg/m^3]");
diff --git a/dumux/freeflow/rans/zeroeq/volumevariables.hh b/dumux/freeflow/rans/zeroeq/volumevariables.hh
index d6555dd72f..fca4665c89 100644
--- a/dumux/freeflow/rans/zeroeq/volumevariables.hh
+++ b/dumux/freeflow/rans/zeroeq/volumevariables.hh
@@ -105,9 +105,12 @@ public:
const Element &element,
const SubControlVolume& scv)
{
+ using std::abs;
using std::exp;
+ using std::sqrt;
Scalar kinematicEddyViscosity = 0.0;
const Scalar karmanConstant = 0.41; // TODO make karman constant a property
+ Scalar velGrad = abs(asImp_().velocityGradients()[0][1]); // TODO: flow and wallnormalaxis
if (eddyViscosityModel_ == Indices::noEddyViscosityModel)
{
// kinematicEddyViscosity = 0.0
@@ -115,14 +118,13 @@ public:
else if (eddyViscosityModel_ == Indices::prandtl)
{
Scalar mixingLength = karmanConstant * asImp_().wallDistance();
- Scalar velGrad = asImp_().velocityGradients()[0][1]; // TODO: flow and wallnormalaxis
kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
}
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
+ * (1.0 - exp(-asImp_().yPlus() / 26.0))
+ / sqrt(1.0 - exp(-0.26 * asImp_().yPlus()));
kinematicEddyViscosity = mixingLength * mixingLength * velGrad;
}
else
diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh
index 8a9ba994d3..8da07e74df 100644
--- a/test/freeflow/rans/pipelauferproblem.hh
+++ b/test/freeflow/rans/pipelauferproblem.hh
@@ -112,7 +112,7 @@ public:
*/
// \{
- const bool isOnWall(GlobalPosition &globalPos) const
+ bool isOnWall(const GlobalPosition &globalPos) const
{
Scalar localEps_ = 1e-6; // cannot use the epsilon, because ParentType is initialized first
return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + localEps_
@@ -180,14 +180,7 @@ public:
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- PrimaryVariables values = initialAtPos(globalPos);
-
- if(isInlet(globalPos))
- {
- values[velocityXIdx] = inletVelocity_;
- }
-
- return values;
+ return initialAtPos(globalPos);
}
// \}
@@ -204,13 +197,12 @@ public:
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
- PrimaryVariables values;
+ PrimaryVariables values(0.0);
values[pressureIdx] = 1.0e+5;
- if(isInlet(globalPos))
+ if(!isOnWall(globalPos))
{
values[velocityXIdx] = inletVelocity_;
}
- values[velocityYIdx] = 0.0;
return values;
}
diff --git a/test/freeflow/rans/test_pipe_laufer.input b/test/freeflow/rans/test_pipe_laufer.input
index 46ac807569..c3e018c3fd 100644
--- a/test/freeflow/rans/test_pipe_laufer.input
+++ b/test/freeflow/rans/test_pipe_laufer.input
@@ -1,5 +1,5 @@
[TimeLoop]
-DtInitial = 1 # [s]
+DtInitial = 1e-1 # [s]
TEnd = 100 # [s]
[Grid]
@@ -23,6 +23,7 @@ EddyViscosityModel = 1
[Newton]
MaxSteps = 10
MaxRelativeShift = 1e-5
+TargetSteps = 8
[Vtk]
AddVelocity = true
--
GitLab