Skip to content

GitLab

dumux
dumux
Commit 8aaabe31 authored by Ned Coltman's avatar Ned Coltman Committed by Melanie Lipp
Browse files
Options

[rans][zeroeq] Clean up roughness lengths and sandgrain roughness

parent 6b15f9a1
Hide whitespace changes
Inline Side-by-side
Showing with 64 additions and 70 deletions
dumux/freeflow/rans/problem.hh
View file @ 8aaabe31
......@@ -114,7 +114,6 @@ public:
flowNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedFlowNormalAxis_);
wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedWallNormalAxis_);
kinematicViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
sandGrainRoughness_.resize(this->gridGeometry().elementMapper().size(), 0.0);
// store the element indicies for all elements with an intersection on the wall
std::vector<unsigned int> wallElementIndicies;
......@@ -183,7 +182,6 @@ public:
wallElementIdx_[elementIdx] = wallElementIndicies[i];
if ( !(hasParam("RANS.WallNormalAxis")) )
wallNormalAxis_[elementIdx] = wallFaceNormalAxis[i];
sandGrainRoughness_[elementIdx] = asImp_().sandGrainRoughnessAtPos(wallPositions[i][0]);
}
}
}
......@@ -487,16 +485,6 @@ public:
return hasChannelGeometry;
}
/*!
* \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;
}
/*!
* \brief Returns the Karman constant
*/
......@@ -598,7 +586,6 @@ private:
std::vector<Scalar> vorticityTensorScalarProduct_;
std::vector<Scalar> kinematicViscosity_;
std::vector<Scalar> sandGrainRoughness_;
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
......
dumux/freeflow/rans/zeroeq/problem.hh
View file @ 8aaabe31
......@@ -105,52 +105,9 @@ public:
{
ParentType::updateDynamicWallProperties(curSol);
// calculate additional roughness
bool printedRangeWarning = false;
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
using std::sqrt;
using std::exp;
const int dofIdx = scv.dofIndex();
// construct a privars object from the cell center solution vector
const auto& cellCenterPriVars = curSol[GridGeometry::cellCenterIdx()][dofIdx];
PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
auto elemSol = elementSolution<typename GridGeometry::LocalView>(std::move(priVars));
VolumeVariables volVars;
volVars.update(elemSol, asImp_(), element, scv);
Scalar ksPlus = this->sandGrainRoughness_[elementIdx] * volVars.uStar() / volVars.kinematicViscosity();
if (ksPlus > 0 && eddyViscosityModel_.compare("baldwinLomax") == 0)
{
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 " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
<< " 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 " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
<< " is not in the valid range (ksPlus > 4.535) and now set to 0.0"<< std::endl;
ksPlus = 0.0;
printedRangeWarning = true;
}
additionalRoughnessLength_[elementIdx] = 0.9 / (volVars.uStar() / volVars.kinematicViscosity())
* (sqrt(ksPlus) - ksPlus * exp(-ksPlus / 6.0));
}
}
// correct roughness lengths if a sand grain roughness is specified
if (hasParam("Problem.SandGrainRoughness"))
calculateRoughnessLength_(curSol);
// update routine for specfic models
if (eddyViscosityModel_.compare("baldwinLomax") == 0)
......@@ -187,7 +144,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
unsigned int flowNormalAxis = this->flowNormalAxis(elementIdx);
unsigned int wallNormalAxis = this->wallNormalAxis(elementIdx);
......@@ -211,7 +168,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
Scalar maxVelocityNorm = 0.0;
Scalar minVelocityNorm = 0.0;
......@@ -238,7 +195,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
// checks if sign switches, by multiplication
Scalar check = kinematicEddyViscosityDifference[asImp_().wallElementIndex(elementIdx)] * kinematicEddyViscosityDifference[elementIdx];
......@@ -253,7 +210,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityInner[elementIdx];
if (effectiveWallDistance >= switchingPosition[asImp_().wallElementIndex(elementIdx)])
......@@ -266,9 +223,64 @@ public:
public:
std::string eddyViscosityModel_;
std::vector<Scalar> kinematicEddyViscosity_;
std::vector<Scalar> additionalRoughnessLength_;
int additionalRoughnessLength(const int elementIdx) const
{ return additionalRoughnessLength_[elementIdx]; }
private:
void calculateRoughnessLength_(const SolutionVector& curSol)
{
bool printedRangeWarning = false;
for (const auto& element : elements(this->gridGeometry().gridView()))
{
static const Scalar sandGrainRoughness = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.SandGrainRoughness", 0.0);
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
using std::sqrt;
using std::exp;
const int dofIdx = scv.dofIndex();
// construct a privars object from the cell center solution vector
const auto& cellCenterPriVars = curSol[GridGeometry::cellCenterIdx()][dofIdx];
PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
auto elemSol = elementSolution<typename GridGeometry::LocalView>(std::move(priVars));
VolumeVariables volVars;
volVars.update(elemSol, asImp_(), element, scv);
Scalar ksPlus = sandGrainRoughness * volVars.uStar() / volVars.kinematicViscosity();
if (ksPlus > 0 && eddyViscosityModel().compare("baldwinLomax") == 0)
{
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 " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
<< " 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 " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
<< " is not in the valid range (ksPlus > 4.535) and now set to 0.0"<< std::endl;
ksPlus = 0.0;
printedRangeWarning = true;
}
additionalRoughnessLength_[elementIdx] = 0.9 / (volVars.uStar() / volVars.kinematicViscosity())
* (sqrt(ksPlus) - ksPlus * exp(-ksPlus / 6.0));
}
}
}
std::vector<Scalar> additionalRoughnessLength_;
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
......
dumux/freeflow/rans/zeroeq/volumevariables.hh
View file @ 8aaabe31
......@@ -81,7 +81,7 @@ public:
const SubControlVolume& scv)
{
RANSParentType::updateRANSProperties(elemSol, problem, element, scv);
additionalRoughnessLength_ = problem.additionalRoughnessLength_[RANSParentType::elementIdx()];
additionalRoughnessLength_ = problem.additionalRoughnessLength(RANSParentType::elementIdx());
yPlusRough_ = wallDistanceRough() * RANSParentType::uStar() / RANSParentType::kinematicViscosity();
RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity(elemSol, problem, element, scv, problem.eddyViscosityModel_));
RANSParentType::calculateEddyDiffusivity(problem);
......
test/freeflow/rans/problem.hh
View file @ 8aaabe31
......@@ -145,7 +145,6 @@ public:
inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
inletTemperature_ = getParam<Scalar>("Problem.InletTemperature", 283.15);
wallTemperature_ = getParam<Scalar>("Problem.WallTemperature", 303.15);
sandGrainRoughness_ = getParam<Scalar>("Problem.SandGrainRoughness", 0.0);
FluidSystem::init();
Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
......@@ -186,9 +185,6 @@ public:
|| globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_;
}
Scalar sandGrainRoughnessAtPos(const GlobalPosition &globalPos) const
{ return sandGrainRoughness_; }
/*!
* \brief Returns the temperature [K] within the domain for the isothermal model.
*/
......@@ -441,7 +437,6 @@ private:
Scalar inletVelocity_;
Scalar inletTemperature_;
Scalar wallTemperature_;
Scalar sandGrainRoughness_;
Scalar initializationTime_;
Scalar viscosityTilde_;
Scalar turbulentKineticEnergy_;
......
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