Commit 6b15f9a1 authored by Ned Coltman's avatar Ned Coltman Committed by Melanie Lipp

[rans][problem] store private variables and output properly

parent 19ad1052
......@@ -40,9 +40,9 @@ struct RANSIOFields
{
NavierStokesIOFields::initOutputModule(out);
static constexpr auto dim = decltype(std::declval<typename OutputModule::VolumeVariables>().velocity())::dimension;
static constexpr auto dim = decltype(std::declval<typename OutputModule::VolumeVariables>().ccVelocityVector())::dimension;
out.addVolumeVariable([](const auto& v){ return v.velocity()[0] / v.velocityMaximum()[0]; }, "v_x/v_x,max");
out.addVolumeVariable([](const auto& v){ return v.ccVelocityVector()[0] / v.profileVelocityMaximum()[0]; }, "v_x/v_x,max");
out.addVolumeVariable([](const auto& v){ return v.velocityGradients()[0]; }, "dv_x/dx_");
if (dim > 1)
out.addVolumeVariable([](const auto& v){ return v.velocityGradients()[1]; }, "dv_y/dx_");
......
......@@ -123,11 +123,13 @@ public:
for (unsigned int dimIdx = 0; dimIdx < Grid::dimension; ++dimIdx)
{
const unsigned int neighborIndex0 = ParentType::neighborIndex(elementIdx, dimIdx, 0);
const unsigned int neighborIndex1 = ParentType::neighborIndex(elementIdx, dimIdx, 1);
// calculate cell-centered turbulentEddyViscosity gradient
storedViscosityTildeGradient_[elementIdx][dimIdx]
= (storedViscosityTilde_[ParentType::neighborIdx_[elementIdx][dimIdx][1]]
- storedViscosityTilde_[ParentType::neighborIdx_[elementIdx][dimIdx][0]])
/ (ParentType::cellCenter_[ParentType::neighborIdx_[elementIdx][dimIdx][1]][dimIdx]
- ParentType::cellCenter_[ParentType::neighborIdx_[elementIdx][dimIdx][0]][dimIdx]);
= (storedViscosityTilde_[neighborIndex1] - storedViscosityTilde_[neighborIndex0])
/ (ParentType::cellCenter(neighborIndex1)[dimIdx] - ParentType::cellCenter(neighborIndex0)[dimIdx]);
}
auto fvGeometry = localView(this->gridGeometry());
......@@ -140,13 +142,13 @@ public:
// face Value
Scalar dirichletViscosityTilde = asImp_().dirichlet(element, scvf)[Indices::viscosityTildeIdx];
unsigned int neighborIdx = ParentType::neighborIdx_[elementIdx][normDim][0];
if (scvf.center()[normDim] < ParentType::cellCenter_[elementIdx][normDim])
neighborIdx = ParentType::neighborIdx_[elementIdx][normDim][1];
unsigned int neighborIndex = ParentType::neighborIndex(elementIdx, normDim, 0);
if (scvf.center()[normDim] < ParentType::cellCenter(elementIdx)[normDim])
neighborIndex = ParentType::neighborIndex(elementIdx, normDim, 1);
storedViscosityTildeGradient_[elementIdx][normDim]
= (storedViscosityTilde_[neighborIdx] - dirichletViscosityTilde)
/ (ParentType::cellCenter_[neighborIdx][normDim] - scvf.center()[normDim]);
= (storedViscosityTilde_[neighborIndex] - dirichletViscosityTilde)
/ (ParentType::cellCenter(neighborIndex)[normDim] - scvf.center()[normDim]);
}
}
}
......
......@@ -87,8 +87,8 @@ public:
viscosityTilde_ = elemSol[0][Indices::viscosityTildeIdx];
storedViscosityTilde_ = problem.storedViscosityTilde_[RANSParentType::elementIdx()];
storedViscosityTildeGradient_ = problem.storedViscosityTildeGradient_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()];
vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct(RANSParentType::elementIdx());
if (problem.useStoredEddyViscosity_)
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]);
else
......@@ -274,8 +274,8 @@ protected:
Scalar viscosityTilde_ = 0.0;
Scalar storedViscosityTilde_ = 0.0;
DimVector storedViscosityTildeGradient_ = DimVector(0.0);
Scalar stressTensorScalarProduct_ = 0.0;
Scalar vorticityTensorScalarProduct_ = 0.0;
Scalar stressTensorScalarProduct_;
Scalar vorticityTensorScalarProduct_;
};
} // end namespace Dumux
......
......@@ -108,7 +108,6 @@ public:
neighborIdx_.resize(this->gridGeometry().elementMapper().size());
cellCenter_.resize(this->gridGeometry().elementMapper().size(), GlobalPosition(0.0));
velocity_.resize(this->gridGeometry().elementMapper().size(), DimVector(0.0));
velocityMaximum_.resize(this->gridGeometry().elementMapper().size(), DimVector(0.0));
velocityGradients_.resize(this->gridGeometry().elementMapper().size(), DimMatrix(0.0));
stressTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
vorticityTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
......@@ -207,16 +206,13 @@ public:
unsigned int neighborIdx = this->gridGeometry().elementMapper().index(intersection.outside());
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
if (abs(cellCenter_[elementIdx][dimIdx] - cellCenter_[neighborIdx][dimIdx]) > 1e-8)
if (abs(cellCenter(elementIdx)[dimIdx] - cellCenter(neighborIdx)[dimIdx]) > 1e-8)
{
if (cellCenter_[elementIdx][dimIdx] > cellCenter_[neighborIdx][dimIdx])
{
if (cellCenter(elementIdx)[dimIdx] > cellCenter(neighborIdx)[dimIdx])
neighborIdx_[elementIdx][dimIdx][0] = neighborIdx;
}
if (cellCenter_[elementIdx][dimIdx] < cellCenter_[neighborIdx][dimIdx])
{
if (cellCenter(elementIdx)[dimIdx] < cellCenter(neighborIdx)[dimIdx])
neighborIdx_[elementIdx][dimIdx][1] = neighborIdx;
}
}
}
}
......@@ -268,40 +264,34 @@ public:
// calculate cell-center-averaged velocity gradients, maximum, and minimum values
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
unsigned int wallElementIdx = wallElementIdx_[elementIdx];
const unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
const unsigned int wallElementIdx = wallElementIdx_[elementIdx];
Scalar maxVelocity = 0.0;
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
{
const unsigned int neighborIndex0 = neighborIndex(elementIdx, dimIdx, 0);
const unsigned int neighborIndex1 = neighborIndex(elementIdx, dimIdx, 1);
velocityGradients_[elementIdx][velIdx][dimIdx]
= (velocity_[neighborIdx_[elementIdx][dimIdx][1]][velIdx]
- velocity_[neighborIdx_[elementIdx][dimIdx][0]][velIdx])
/ (cellCenter_[neighborIdx_[elementIdx][dimIdx][1]][dimIdx]
- cellCenter_[neighborIdx_[elementIdx][dimIdx][0]][dimIdx]);
if (abs(cellCenter_[neighborIdx_[elementIdx][dimIdx][1]][dimIdx]
- cellCenter_[neighborIdx_[elementIdx][dimIdx][0]][dimIdx]) < 1e-8)
= (ccVelocity(neighborIndex1, velIdx) - ccVelocity(neighborIndex0, velIdx))
/ (cellCenter(neighborIndex1)[dimIdx] - cellCenter(neighborIndex0)[dimIdx]);
if (abs(cellCenter(neighborIndex1)[dimIdx] - cellCenter(neighborIndex0)[dimIdx]) < 1e-8)
velocityGradients_[elementIdx][velIdx][dimIdx] = 0.0;
}
if (abs(velocity_[elementIdx][dimIdx]) > abs(velocityMaximum_[wallElementIdx][dimIdx]))
{
velocityMaximum_[wallElementIdx][dimIdx] = velocity_[elementIdx][dimIdx];
}
if (abs(velocity_[elementIdx][dimIdx]) < abs(velocityMinimum_[wallElementIdx][dimIdx]))
{
velocityMinimum_[wallElementIdx][dimIdx] = velocity_[elementIdx][dimIdx];
}
if (abs(ccVelocity(elementIdx, dimIdx)) > abs(velocityMaximum_[wallElementIdx][dimIdx]))
velocityMaximum_[wallElementIdx][dimIdx] = ccVelocity(elementIdx, dimIdx);
if (0 <= flowNormalAxis && flowNormalAxis < dim)
{
flowNormalAxis_[elementIdx] = flowNormalAxis;
}
else if (abs(maxVelocity) < abs(velocity_[elementIdx][dimIdx]))
if (abs(ccVelocity(elementIdx, dimIdx)) < abs(velocityMinimum_[wallElementIdx][dimIdx]))
velocityMinimum_[wallElementIdx][dimIdx] = ccVelocity(elementIdx, dimIdx);
if ((hasParam("RANS.FlowNormalAxis") != 1) && (maxVelocity) < abs(ccVelocity(elementIdx, dimIdx)))
{
maxVelocity = abs(velocity_[elementIdx][dimIdx]);
maxVelocity = abs(ccVelocity(elementIdx, dimIdx));
flowNormalAxis_[elementIdx] = dimIdx;
}
}
......@@ -321,13 +311,13 @@ public:
Scalar dirichletVelocity = asImp_().dirichlet(element, scvf)[Indices::velocity(velIdx)];
unsigned int neighborIdx = neighborIdx_[elementIdx][scvfNormDim][0];
if (scvf.center()[scvfNormDim] < cellCenter_[elementIdx][scvfNormDim])
neighborIdx = neighborIdx_[elementIdx][scvfNormDim][1];
unsigned int neighborIdx = neighborIndex(elementIdx, scvfNormDim, 0);
if (scvf.center()[scvfNormDim] < cellCenter(elementIdx)[scvfNormDim])
neighborIdx = neighborIndex(elementIdx, scvfNormDim, 1);
velocityGradients_[elementIdx][velIdx][scvfNormDim]
= (velocity_[neighborIdx][velIdx] - dirichletVelocity)
/ (cellCenter_[neighborIdx][scvfNormDim] - scvf.center()[scvfNormDim]);
= (ccVelocity(neighborIdx, velIdx) - dirichletVelocity)
/ (cellCenter(neighborIdx)[scvfNormDim] - scvf.center()[scvfNormDim]);
}
}
......@@ -346,12 +336,12 @@ public:
unsigned int normalNormDim = lateralFace.directionIndex();
if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isBeaversJoseph(Indices::velocity(velIdx))))
{
unsigned int neighborIdx = neighborIdx_[elementIdx][normalNormDim][0];
if (lateralFace.center()[normalNormDim] < cellCenter_[elementIdx][normalNormDim])
neighborIdx = neighborIdx_[elementIdx][normalNormDim][1];
unsigned int neighborIdx = neighborIndex(elementIdx, normalNormDim, 0);
if (lateralFace.center()[normalNormDim] < cellCenter(elementIdx)[normalNormDim])
neighborIdx = neighborIndex(elementIdx, normalNormDim, 1);
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
bjsVelocityAverage[normalNormDim] += ParentType::beaversJosephVelocity(element, scv, scvf, lateralFace, velocity_[elementIdx][velIdx], 0.0);
bjsVelocityAverage[normalNormDim] += ParentType::beaversJosephVelocity(element, scv, scvf, lateralFace, ccVelocity(elementIdx, velIdx), 0.0);
if (bjsNumFaces[normalNormDim] > 0 && neighborIdx != bjsNeighbor[normalNormDim])
DUNE_THROW(Dune::InvalidStateException, "Two different neighborIdx should not occur");
bjsNeighbor[normalNormDim] = neighborIdx;
......@@ -368,8 +358,8 @@ public:
bjsVelocityAverage[dirIdx] /= bjsNumFaces[dirIdx];
velocityGradients_[elementIdx][velIdx][dirIdx]
= (velocity_[neighborIdx][velIdx] - bjsVelocityAverage[dirIdx])
/ (cellCenter_[neighborIdx][dirIdx] - normalNormCoordinate[dirIdx]);
= (ccVelocity(neighborIdx, velIdx) - bjsVelocityAverage[dirIdx])
/ (cellCenter(neighborIdx)[dirIdx] - normalNormCoordinate[dirIdx]);
}
}
......@@ -385,8 +375,8 @@ public:
{
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
{
stressTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementIdx][dimIdx][velIdx]
+ 0.5 * velocityGradients_[elementIdx][velIdx][dimIdx];
stressTensor[dimIdx][velIdx] = 0.5 * velocityGradient(elementIdx, dimIdx, velIdx)
+ 0.5 * velocityGradient(elementIdx, velIdx, dimIdx);
}
}
stressTensorScalarProduct_[elementIdx] = 0.0;
......@@ -403,8 +393,8 @@ public:
{
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
{
vorticityTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementIdx][dimIdx][velIdx]
- 0.5 * velocityGradients_[elementIdx][velIdx][dimIdx];
vorticityTensor[dimIdx][velIdx] = 0.5 * velocityGradient(elementIdx, dimIdx, velIdx)
- 0.5 * velocityGradient(elementIdx, velIdx, dimIdx);
}
}
vorticityTensorScalarProduct_[elementIdx] = 0.0;
......@@ -541,7 +531,6 @@ public:
return turbulentSchmidtNumber;
}
public:
int wallNormalAxis(const int elementIdx) const
{ return wallNormalAxis_[elementIdx]; }
......@@ -554,9 +543,42 @@ public:
Scalar wallDistance(const int elementIdx) const
{ return wallDistance_[elementIdx]; }
GlobalPosition cellCenter(const int elementIdx) const
{ return cellCenter_[elementIdx]; }
unsigned int neighborIndex(const int elementIdx, const int dimIdx, const int sideIdx) const
{ return neighborIdx_[elementIdx][dimIdx][sideIdx];}
DimVector ccVelocityVector(const int elementIdx) const
{ return velocity_[elementIdx]; }
Scalar ccVelocity(const int elementIdx, const int dimIdx) const
{ return velocity_[elementIdx][dimIdx]; }
DimVector profileVelocityMaximum(const int elementIdx) const
{ return velocityMaximum_[elementIdx]; }
DimVector profileVelocityMinimum(const int elementIdx) const
{ return velocityMinimum_[elementIdx]; }
DimMatrix velocityGradientTensor(const int elementIdx) const
{ return velocityGradients_[elementIdx]; }
Scalar velocityGradient(const int elementIdx, const int dimIdx, const int velIdx) const
{ return velocityGradients_[elementIdx][dimIdx][velIdx]; }
Scalar stressTensorScalarProduct(const int elementIdx) const
{ return stressTensorScalarProduct_[elementIdx]; }
Scalar vorticityTensorScalarProduct(const int elementIdx) const
{ return vorticityTensorScalarProduct_[elementIdx]; }
Scalar kinematicViscosity(const int elementIdx) const
{ return kinematicViscosity_[elementIdx]; }
bool calledUpdateStaticWallProperties = false;
private:
const int fixedFlowNormalAxis_ = getParam<int>("RANS.FlowNormalAxis", 0);
const int fixedWallNormalAxis_ = getParam<int>("RANS.WallNormalAxis", 1);
......@@ -564,8 +586,9 @@ public:
std::vector<unsigned int> flowNormalAxis_;
std::vector<Scalar> wallDistance_;
std::vector<unsigned int> wallElementIdx_;
std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIdx_;
std::vector<GlobalPosition> cellCenter_;
std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIdx_;
std::vector<DimVector> velocity_;
std::vector<DimVector> velocityMaximum_;
std::vector<DimVector> velocityMinimum_;
......@@ -577,7 +600,6 @@ public:
std::vector<Scalar> kinematicViscosity_;
std::vector<Scalar> sandGrainRoughness_;
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
......
......@@ -204,7 +204,7 @@ public:
unsigned int wallElementIdx = asImp_().wallElementIndex(elementIdx);
unsigned int matchingPointIdx = matchingPointIdx_[wallElementIdx];
return wallElementIdx == matchingPointIdx ? yPlusNominal(elementIdx) < yPlusThreshold_
: yPlus(elementIdx) < yPlusThreshold_;
: yPlus(elementIdx) < yPlusThreshold_;
}
/*!
......
......@@ -87,11 +87,11 @@ public:
dissipation_ = elemSol[0][Indices::dissipationIdx];
storedDissipation_ = problem.storedDissipation_[RANSParentType::elementIdx()];
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
const Scalar uStarNominal = problem.uStarNominal(RANSParentType::elementIdx());
const auto flowNormalAxis = problem.flowNormalAxis_[RANSParentType::elementIdx()];
yPlusNominal_ = RANSParentType::wallDistance() * uStarNominal / problem.kinematicViscosity_[RANSParentType::elementIdx()];
uPlusNominal_ = RANSParentType::velocity()[flowNormalAxis] / uStarNominal;
const auto flowNormalAxis = problem.flowNormalAxis(RANSParentType::elementIdx());
yPlusNominal_ = RANSParentType::wallDistance() * uStarNominal / problem.kinematicViscosity(RANSParentType::elementIdx());
uPlusNominal_ = RANSParentType::ccVelocityVector()[flowNormalAxis] / uStarNominal;
cMu_ = problem.cMu();
if (problem.useStoredEddyViscosity_)
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]);
......
......@@ -121,18 +121,17 @@ public:
for (unsigned int dimIdx = 0; dimIdx < DimVector::dimension; ++dimIdx)
{
unsigned backwardNeighbor = ParentType::neighborIdx_[elementIdx][dimIdx][0];
unsigned forwardNeighbor = ParentType::neighborIdx_[elementIdx][dimIdx][1];
const unsigned neighborIdx0 = ParentType::neighborIndex(elementIdx, dimIdx, 0);
const unsigned neighborIdx1 = ParentType::neighborIndex(elementIdx, dimIdx, 1);
// Cell centered TKE Gradient
storedTurbulentKineticEnergyGradient_[elementIdx][dimIdx]
= (storedTurbulentKineticEnergy_[forwardNeighbor]
- storedTurbulentKineticEnergy_[backwardNeighbor])
/ (ParentType::cellCenter_[forwardNeighbor][dimIdx]
- ParentType::cellCenter_[backwardNeighbor][dimIdx]);
= (storedTurbulentKineticEnergy_[neighborIdx1] - storedTurbulentKineticEnergy_[neighborIdx0])
/ (ParentType::cellCenter(neighborIdx1)[dimIdx] - ParentType::cellCenter(neighborIdx0)[dimIdx]);
// Cell centered Omega Gradient
storedDissipationGradient_[elementIdx][dimIdx]
= (storedDissipation_[forwardNeighbor]
- storedDissipation_[backwardNeighbor])
/ (ParentType::cellCenter_[forwardNeighbor][dimIdx]
- ParentType::cellCenter_[backwardNeighbor][dimIdx]);
= (storedDissipation_[neighborIdx1] - storedDissipation_[neighborIdx0])
/ (ParentType::cellCenter(neighborIdx1)[dimIdx] - ParentType::cellCenter(neighborIdx0)[dimIdx]);
}
}
}
......
......@@ -90,7 +90,7 @@ public:
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementIdx()];
storedDissipationGradient_ = problem.storedDissipationGradient_[RANSParentType::elementIdx()];
storedTurbulentKineticEnergyGradient_ = problem.storedTurbulentKineticEnergyGradient_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
if (problem.useStoredEddyViscosity_)
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]);
else
......
......@@ -86,7 +86,7 @@ public:
dissipationTilde_ = elemSol[0][Indices::dissipationIdx];
storedDissipationTilde_ = problem.storedDissipationTilde_[RANSParentType::elementIdx()];
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
if (problem.useStoredEddyViscosity_)
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]);
else
......
......@@ -94,16 +94,17 @@ public:
// calculate characteristic properties of the turbulent flow
elementIdx_ = problem.gridGeometry().elementMapper().index(element);
wallDistance_ = problem.wallDistance(elementIdx_);
velocity_ = problem.velocity_[elementIdx_];
velocityMaximum_ = problem.velocityMaximum_[wallElementIdx_];
velocityGradients_ = problem.velocityGradients_[elementIdx_];
uStar_ = sqrt(problem.kinematicViscosity_[wallElementIdx_]
* abs(problem.velocityGradients_[wallElementIdx_][flowNormalAxis][wallNormalAxis]));
ccVelocityVector_ = problem.ccVelocityVector(elementIdx_);
profileVelocityMaximum_ = problem.profileVelocityMaximum(problem.wallElementIndex(elementIdx_));
profileVelocityMinimum_ = problem.profileVelocityMinimum(problem.wallElementIndex(elementIdx_));
velocityGradientTensor_ = problem.velocityGradientTensor(elementIdx_);
const auto flowNormalAxis = problem.flowNormalAxis(elementIdx_);
const auto wallNormalAxis = problem.wallNormalAxis(elementIdx_);
uStar_ = sqrt(problem.kinematicViscosity(problem.wallElementIndex(elementIdx_))
* abs(problem.velocityGradient(problem.wallElementIndex(elementIdx_), flowNormalAxis, wallNormalAxis)));
uStar_ = max(uStar_, 1e-10); // zero values lead to numerical problems in some turbulence models
yPlus_ = wallDistance_ * uStar_ / problem.kinematicViscosity_[elementIdx_];
uPlus_ = velocity_[flowNormalAxis] / uStar_;
yPlus_ = wallDistance_ * uStar_ / problem.kinematicViscosity(elementIdx_);
uPlus_ = problem.ccVelocity(elementIdx_, flowNormalAxis) / uStar_;
karmanConstant_ = problem.karmanConstant();
}
......@@ -116,20 +117,26 @@ public:
/*!
* \brief Return the velocity vector \f$\mathrm{[m/s]}\f$ at the control volume center.
*/
DimVector velocity() const
{ return velocity_; }
DimVector ccVelocityVector() const
{ return ccVelocityVector_; }
/*!
* \brief Return the maximum velocity vector \f$\mathrm{[m/s]}\f$ of the wall segment.
*/
DimVector velocityMaximum() const
{ return velocityMaximum_; }
DimVector profileVelocityMaximum() const
{ return profileVelocityMaximum_; }
/*!
* \brief Return the minimum velocity vector \f$\mathrm{[m/s]}\f$ of the wall segment.
*/
DimVector profileVelocityMinimum() const
{ return profileVelocityMinimum_; }
/*!
* \brief Return the velocity gradients \f$\mathrm{[1/s]}\f$ at the control volume center.
*/
DimMatrix velocityGradients() const
{ return velocityGradients_; }
{ return velocityGradientTensor_; }
/*!
* \brief Return the wall distance \f$\mathrm{[m]}\f$ of the control volume.
......@@ -253,9 +260,10 @@ protected:
Scalar setDynamicEddyViscosity_(Scalar value)
{ return dynamicEddyViscosity_ = value; }
DimVector velocity_;
DimVector velocityMaximum_;
DimMatrix velocityGradients_;
DimVector ccVelocityVector_;
DimVector profileVelocityMaximum_;
DimVector profileVelocityMinimum_;
DimMatrix velocityGradientTensor_;
std::size_t elementIdx_;
Scalar wallDistance_;
Scalar karmanConstant_;
......
......@@ -26,6 +26,7 @@
#include <string>
#include <dune/common/power.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/staggeredfvproblem.hh>
#include <dumux/discretization/localview.hh>
......@@ -190,20 +191,19 @@ public:
unsigned int flowNormalAxis = this->flowNormalAxis(elementIdx);
unsigned int wallNormalAxis = this->wallNormalAxis(elementIdx);
Scalar omegaAbs = abs(this->velocityGradients_[elementIdx][flowNormalAxis][wallNormalAxis]
- this->velocityGradients_[elementIdx][wallNormalAxis][flowNormalAxis]);
Scalar uStar = sqrt(this->kinematicViscosity_[wallElementIdx]
* abs(this->velocityGradients_[wallElementIdx][flowNormalAxis][wallNormalAxis]));
Scalar yPlus = wallDistance * uStar / this->kinematicViscosity_[elementIdx];
Scalar mixingLength = this->karmanConstant() * wallDistance * (1.0 - exp(-yPlus / aPlus));
Scalar omegaAbs = abs(this->velocityGradient(elementIdx, flowNormalAxis, wallNormalAxis)
- this->velocityGradient(elementIdx, wallNormalAxis, flowNormalAxis));
Scalar uStar = sqrt(this->kinematicViscosity(asImp_().wallElementIndex(elementIdx))
* abs(this->velocityGradient(asImp_().wallElementIndex(elementIdx), flowNormalAxis, wallNormalAxis)));
Scalar yPlus = effectiveWallDistance * uStar / this->kinematicViscosity(elementIdx);
Scalar mixingLength = this->karmanConstant() * effectiveWallDistance * (1.0 - exp(-yPlus / aPlus));
kinematicEddyViscosityInner[elementIdx] = mixingLength * mixingLength * omegaAbs;
Scalar f = wallDistance * omegaAbs * (1.0 - exp(-yPlus / aPlus));
if (f > storedFMax[wallElementIdx])
Scalar f = effectiveWallDistance * omegaAbs * (1.0 - exp(-yPlus / aPlus));
if (f > storedFMax[asImp_().wallElementIndex(elementIdx)])
{
storedFMax[wallElementIdx] = f;
storedYFMax[wallElementIdx] = wallDistance;
storedFMax[asImp_().wallElementIndex(elementIdx)] = f;
storedYFMax[asImp_().wallElementIndex(elementIdx)] = effectiveWallDistance;
}
}
......@@ -217,16 +217,17 @@ public:
Scalar minVelocityNorm = 0.0;
for (unsigned dimIdx = 0; dimIdx < dim; ++dimIdx)
{
maxVelocityNorm += this->velocityMaximum_[wallElementIdx][dimIdx]
* this->velocityMaximum_[wallElementIdx][dimIdx];
minVelocityNorm += this->velocityMinimum_[wallElementIdx][dimIdx]
* this->velocityMinimum_[wallElementIdx][dimIdx];
maxVelocityNorm += asImp_().profileVelocityMaximum(asImp_().wallElementIndex(elementIdx))[dimIdx]
* asImp_().profileVelocityMaximum(asImp_().wallElementIndex(elementIdx))[dimIdx];
minVelocityNorm += asImp_().profileVelocityMinimum(asImp_().wallElementIndex(elementIdx))[dimIdx]
* asImp_().profileVelocityMinimum(asImp_().wallElementIndex(elementIdx))[dimIdx];
}
Scalar deltaU = sqrt(maxVelocityNorm) - sqrt(minVelocityNorm);
Scalar yFMax = storedYFMax[wallElementIdx];
Scalar fMax = storedFMax[wallElementIdx];
Scalar yFMax = storedYFMax[asImp_().wallElementIndex(elementIdx)];
Scalar fMax = storedFMax[asImp_().wallElementIndex(elementIdx)];
Scalar fWake = min(yFMax * fMax, cWake * yFMax * deltaU * deltaU / fMax);
Scalar fKleb = 1.0 / (1.0 + 5.5 * pow(cKleb * wallDistance / yFMax, 6.0));
Scalar fKleb = 1.0 / (1.0 + 5.5 * power(cKleb * effectiveWallDistance / yFMax, 6));
kinematicEddyViscosityOuter[elementIdx] = k * cCP * fWake * fKleb;
kinematicEddyViscosityDifference[elementIdx]
......@@ -240,11 +241,11 @@ public:
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
// checks if sign switches, by multiplication
Scalar check = kinematicEddyViscosityDifference[wallElementIdx] * kinematicEddyViscosityDifference[elementIdx];
Scalar check = kinematicEddyViscosityDifference[asImp_().wallElementIndex(elementIdx)] * kinematicEddyViscosityDifference[elementIdx];
if (check < 0 // means sign has switched
&& switchingPosition[wallElementIdx] > wallDistance)
&& switchingPosition[asImp_().wallElementIndex(elementIdx)] > effectiveWallDistance)
{
switchingPosition[wallElementIdx] = wallDistance;
switchingPosition[asImp_().wallElementIndex(elementIdx)] = effectiveWallDistance;
}
}
......
......@@ -431,7 +431,7 @@ private:
// For the komega model we set a fixed value for the dissipation
const auto wallDistance = ParentType::wallDistance(elementIdx);
using std::pow;
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx]
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
/ (ParentType::betaOmega() * wallDistance * wallDistance);
return values;
}
......
......@@ -445,9 +445,9 @@ private:
{
static_assert(ModelTraits::turbulenceModel() == TurbulenceModel::komega, "Only valid for Komega");
// For the komega model we set a fixed value for the dissipation
const auto wallDistance = ParentType::wallDistance_[elementIdx];
const auto wallDistance = ParentType::wallDistance(elementIdx);
using std::pow;
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx]
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity(elementIdx)
/ (ParentType::betaOmega() * wallDistance * wallDistance);
return values;
}
......
Markdown is supported
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