Skip to content

GitLab

Commit c0338acb authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[rans][zeroeq] First version of velocity gradients

parent 5470609d
Hide whitespace changes
Inline Side-by-side
dumux/freeflow/rans/problem.hh
View file @ c0338acb
......@@ -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_;
......
dumux/freeflow/rans/volumevariables.hh
View file @ c0338acb
......@@ -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);
......
dumux/freeflow/rans/vtkoutputfields.hh
View file @ c0338acb
......@@ -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]");
......
dumux/freeflow/rans/zeroeq/volumevariables.hh
View file @ c0338acb
......@@ -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
......
test/freeflow/rans/pipelauferproblem.hh
View file @ c0338acb
......@@ -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;
}
......
test/freeflow/rans/test_pipe_laufer.input
View file @ c0338acb
[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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment