Commit 16cd0d39 authored by Melanie Lipp's avatar Melanie Lipp

Merge branch 'feature/fix_walldistance' into 'master'

Feature/fix walldistance

See merge request !2224
parents fb2a953c 6a01e63e
......@@ -333,6 +333,9 @@ private:
// parameters in the mpfa group
defaultParams["MPFA.Q"] = "0.0";
// parameters in the RANS group
defaultParams["RANS.IsFlatWallBounded"] = "false";
// merge the global default tree but do not overwrite if the parameter already exists
mergeTree_(params, defaultParams, false);
}
......
......@@ -142,6 +142,8 @@ class FreeFlowStaggeredSubControlVolumeFace
public:
using GlobalPosition = typename T::GlobalPosition;
static constexpr int numCornersPerFace = 2 * (dim - 1);
//! State the traits public and thus export all types
using Traits = T;
......
......@@ -188,9 +188,7 @@ public:
* \param compIdx the index of the component
*/
Scalar massFraction(int phaseIdx, int compIdx) const
{
return fluidState_.massFraction(0, compIdx);
}
{ return fluidState_.massFraction(0, compIdx); }
/*!
* \brief Returns the mass fraction of a component in the phase \f$\mathrm{[-]}\f$
......@@ -198,9 +196,7 @@ public:
* \param compIdx the index of the component
*/
Scalar massFraction(int compIdx) const
{
return fluidState_.massFraction(0, compIdx);
}
{ return fluidState_.massFraction(0, compIdx); }
/*!
* \brief Returns the mole fraction of a component in the phase \f$\mathrm{[-]}\f$
......@@ -209,9 +205,7 @@ public:
* \param compIdx the index of the component
*/
Scalar moleFraction(int phaseIdx, int compIdx) const
{
return fluidState_.moleFraction(0, compIdx);
}
{ return fluidState_.moleFraction(0, compIdx); }
/*!
* \brief Returns the mole fraction of a component in the phase \f$\mathrm{[-]}\f$
......@@ -219,17 +213,13 @@ public:
* \param compIdx the index of the component
*/
Scalar moleFraction(int compIdx) const
{
return fluidState_.moleFraction(0, compIdx);
}
{ return fluidState_.moleFraction(0, compIdx); }
/*!
* \brief Returns the mass density of a given phase \f$\mathrm{[kg/m^3]}\f$
*/
Scalar molarDensity(int phaseIdx = 0) const
{
return fluidState_.molarDensity(0);
}
{ return fluidState_.molarDensity(0); }
/*!
* \brief Returns the molar mass of a given component.
......@@ -237,9 +227,7 @@ public:
* \param compIdx the index of the component
*/
Scalar molarMass(int compIdx) const
{
return FluidSystem::molarMass(compIdx);
}
{ return FluidSystem::molarMass(compIdx); }
/*!
* \brief Returns the average molar mass \f$\mathrm{[kg/mol]}\f$ of the fluid phase.
......
......@@ -40,9 +40,11 @@ struct RANSIOFields
{
NavierStokesIOFields::initOutputModule(out);
static constexpr auto dim = decltype(std::declval<typename OutputModule::VolumeVariables>().velocity())::dimension;
static const bool isFlatWallBounded = getParamFromGroup<bool>(out.paramGroup(), "RANS.IsFlatWallBounded");
out.addVolumeVariable([](const auto& v){ return v.velocity()[0] / v.velocityMaximum()[0]; }, "v_x/v_x,max");
static constexpr auto dim = decltype(std::declval<typename OutputModule::VolumeVariables>().ccVelocityVector())::dimension;
out.addVolumeVariable([](const auto& v){ return v.ccVelocityVector()[0] / v.velocityMaximum()[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_");
......@@ -52,8 +54,11 @@ struct RANSIOFields
out.addVolumeVariable([](const auto& v){ return v.viscosity() / v.density(); }, "nu");
out.addVolumeVariable([](const auto& v){ return v.kinematicEddyViscosity(); }, "nu_t");
out.addVolumeVariable([](const auto& v){ return v.wallDistance(); }, "l_w");
out.addVolumeVariable([](const auto& v){ return v.yPlus(); }, "y^+");
out.addVolumeVariable([](const auto& v){ return v.uPlus(); }, "u^+");
if (isFlatWallBounded)
{
out.addVolumeVariable([](const auto& v){ return v.yPlus(); }, "y^+");
out.addVolumeVariable([](const auto& v){ return v.uPlus(); }, "u^+");
}
}
//! return the names of the primary variables
......
......@@ -68,10 +68,7 @@ public:
//! The constructor sets the gravity, if desired by the user.
RANSProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
: ParentType(gridGeometry, paramGroup)
{
useStoredEddyViscosity_ = getParamFromGroup<bool>(this->paramGroup(),
"RANS.UseStoredEddyViscosity", false);
}
{ }
/*!
* \brief Correct size of the static (solution independent) wall variables
......@@ -116,49 +113,65 @@ public:
}
}
// calculate cell-center-averaged velocity gradients, maximum, and minimum values
// calculate cell-center-averaged gradient
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
const unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
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 (viscosityTilde) 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]);
}
// Adjust for dirichlet boundary conditions
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scvf : scvfs(fvGeometry))
{
unsigned int normDim = scvf.directionIndex();
const unsigned int normDim = scvf.directionIndex();
if (scvf.boundary() && asImp_().boundaryTypes(element, scvf).isDirichlet(Indices::viscosityTildeIdx))
{
// 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]);
}
}
}
}
public:
bool useStoredEddyViscosity() const
{
static const bool useStoredEddyViscosity = getParamFromGroup<bool>(this->paramGroup(), "RANS.UseStoredEddyViscosity", false);
return useStoredEddyViscosity;
}
Scalar storedDynamicEddyViscosity(const int elementIdx) const
{ return storedDynamicEddyViscosity_[elementIdx]; }
Scalar storedViscosityTilde(const int elementIdx) const
{ return storedViscosityTilde_[elementIdx]; }
DimVector storedViscosityTildeGradient(const int elementIdx) const
{ return storedViscosityTildeGradient_[elementIdx]; }
private:
std::vector<Scalar> storedDynamicEddyViscosity_;
std::vector<Scalar> storedViscosityTilde_;
std::vector<DimVector> storedViscosityTildeGradient_;
bool useStoredEddyViscosity_;
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
......
......@@ -85,12 +85,12 @@ public:
{
RANSParentType::updateRANSProperties(elemSol, problem, element, scv);
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()];
if (problem.useStoredEddyViscosity_)
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]);
storedViscosityTilde_ = problem.storedViscosityTilde(RANSParentType::elementIdx());
storedViscosityTildeGradient_ = problem.storedViscosityTildeGradient(RANSParentType::elementIdx());
stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct(RANSParentType::elementIdx());
if (problem.useStoredEddyViscosity())
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity(RANSParentType::elementIdx()));
else
RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity());
RANSParentType::calculateEddyDiffusivity(problem);
......@@ -101,41 +101,31 @@ public:
* \brief Returns the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$.
*/
Scalar calculateEddyViscosity()
{
return viscosityTilde() * fv1() * RANSParentType::density();
}
{ return viscosityTilde() * fv1() * RANSParentType::density(); }
/*!
* \brief Returns the viscosity parameter \f$ m^2/s \f$
*/
Scalar viscosityTilde() const
{
return viscosityTilde_;
}
{ return viscosityTilde_; }
/*!
* \brief Returns the viscosity parameter from the last iteration \f$ m^2/s \f$
*/
Scalar storedViscosityTilde() const
{
return storedViscosityTilde_;
}
{ return storedViscosityTilde_; }
/*!
* \brief Returns the gradient of the viscosity parameter
*/
DimVector storedViscosityTildeGradient() const
{
return storedViscosityTildeGradient_;
}
{ return storedViscosityTildeGradient_; }
/*!
* \brief Returns the scalar product of the stress tensor
*/
Scalar stressTensorScalarProduct() const
{
return stressTensorScalarProduct_;
}
{ return stressTensorScalarProduct_; }
/*!
* \brief Returns damping function for the eddy viscosity
......@@ -274,8 +264,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
......
This diff is collapsed.
......@@ -85,18 +85,21 @@ public:
//! The constructor sets the gravity, if desired by the user.
RANSProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
: ParentType(gridGeometry, paramGroup)
{
yPlusThreshold_ = getParamFromGroup<Scalar>(this->paramGroup(), "KEpsilon.YPlusThreshold", 30);
useStoredEddyViscosity_ = getParamFromGroup<bool>(this->paramGroup(), "RANS.UseStoredEddyViscosity", false);
}
{ }
/*!
* \brief Correct size of the static (solution independent) wall variables
*/
void updateStaticWallProperties()
{
ParentType::updateStaticWallProperties();
if (!ParentType::isFlatWallBounded())
{
DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, k-epsilon models should only be used for "
<< " wall bounded flows with flat channel geometries. "
<< "\n If your geometry is a flat channel, please set the runtime parameter RANS.IsFlatWallBounded to true. \n");
}
ParentType::updateStaticWallProperties();
// update size and initial values of the global vectors
matchingPointIdx_.resize(this->gridGeometry().elementMapper().size(), 0);
storedDensity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
......@@ -144,17 +147,17 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementIdx];
unsigned int neighborIdx0 = asImp_().neighborIdx_[elementIdx][wallNormalAxis][0];
unsigned int neighborIdx1 = asImp_().neighborIdx_[elementIdx][wallNormalAxis][1];
unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
unsigned int neighborIndex0 = asImp_().neighborIndex(elementIdx, wallNormalAxis, 0);
unsigned int neighborIndex1 = asImp_().neighborIndex(elementIdx, wallNormalAxis, 1);
numElementsInNearWallRegion = inNearWallRegion(elementIdx)
? numElementsInNearWallRegion + 1
: numElementsInNearWallRegion + 0;
if ((!inNearWallRegion(elementIdx) && (inNearWallRegion(neighborIdx0) || inNearWallRegion(neighborIdx1)))
|| (!inNearWallRegion(elementIdx) && elementIdx == asImp_().wallElementIdx_[elementIdx])
|| (inNearWallRegion(elementIdx) && (asImp_().wallElementIdx_[neighborIdx0] != asImp_().wallElementIdx_[neighborIdx1])))
if ((!inNearWallRegion(elementIdx) && (inNearWallRegion(neighborIndex0) || inNearWallRegion(neighborIndex1)))
|| (!inNearWallRegion(elementIdx) && elementIdx == asImp_().wallElementIndex(elementIdx))
|| (inNearWallRegion(elementIdx) && (asImp_().wallElementIndex(neighborIndex0) != asImp_().wallElementIndex(neighborIndex1))))
{
matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]] = elementIdx;
matchingPointIdx_[asImp_().wallElementIndex(elementIdx)] = elementIdx;
}
}
std::cout << "numElementsInNearWallRegion: " << numElementsInNearWallRegion << std::endl;
......@@ -174,10 +177,10 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]];
unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx));
Scalar scalingFactor = storedDynamicEddyViscosity_[matchingPointIdx]
/ zeroEqDynamicEddyViscosity_[matchingPointIdx];
Scalar scalingFactor = storedDynamicEddyViscosity(matchingPointIndex)
/ zeroEqDynamicEddyViscosity_[matchingPointIndex];
if (!isMatchingPoint(elementIdx)
&& !std::isnan(scalingFactor) && !std::isinf(scalingFactor))
{
......@@ -187,10 +190,10 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]];
unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx));
if (isMatchingPoint(elementIdx))
{
zeroEqDynamicEddyViscosity_[matchingPointIdx] = storedDynamicEddyViscosity_[matchingPointIdx];
zeroEqDynamicEddyViscosity_[matchingPointIndex] = storedDynamicEddyViscosity(matchingPointIndex);
}
}
}
......@@ -201,33 +204,33 @@ public:
*/
bool inNearWallRegion(unsigned int elementIdx) const
{
unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx];
unsigned int matchingPointIdx = matchingPointIdx_[wallElementIdx];
return wallElementIdx == matchingPointIdx ? yPlusNominal(elementIdx) < yPlusThreshold_
: yPlus(elementIdx) < yPlusThreshold_;
unsigned int wallElementIdx = asImp_().wallElementIndex(elementIdx);
unsigned int matchingPointIndex = matchingPointIdx(wallElementIdx);
return (wallElementIdx == matchingPointIndex) ? yPlusNominal(elementIdx) < yPlusThreshold()
: yPlus(elementIdx) < yPlusThreshold();
}
/*!
* \brief Returns if an element is the matching point
*/
bool isMatchingPoint(unsigned int elementIdx) const
{ return matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]] == elementIdx; }
{ return matchingPointIdx(asImp_().wallElementIndex(elementIdx)) == elementIdx; }
/*!
* \brief Returns the \f$ y^+ \f$ value at an element center
*/
const Scalar yPlus(unsigned int elementIdx) const
{
return asImp_().wallDistance_[elementIdx] * uStar(elementIdx)
/ asImp_().kinematicViscosity_[elementIdx];
return asImp_().wallDistance(elementIdx) * uStar(elementIdx)
/ asImp_().kinematicViscosity(elementIdx);
}
/*!
* \brief Returns the nominal \f$ y^+ \f$ value at an element center
*/
const Scalar yPlusNominal(unsigned int elementIdx) const
{
return asImp_().wallDistance_[elementIdx] * uStarNominal(elementIdx)
/ asImp_().kinematicViscosity_[elementIdx];
return asImp_().wallDistance(elementIdx) * uStarNominal(elementIdx)
/ asImp_().kinematicViscosity(elementIdx);
}
/*!
......@@ -244,15 +247,15 @@ public:
Scalar mixingLength = 0.0;
if (yPlusValue > 0.0)
{
mixingLength = asImp_().karmanConstant() * asImp_().wallDistance_[elementIdx]
mixingLength = asImp_().karmanConstant() * asImp_().wallDistance(elementIdx)
* (1.0 - exp(-yPlusValue / 26.0 ))
/ sqrt(1.0 - exp(-0.26 * yPlusValue));
}
unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementIdx];
unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementIdx];
Scalar velocityGradient = asImp_().velocityGradients_[elementIdx][flowNormalAxis][wallNormalAxis];
return mixingLength * mixingLength * abs(velocityGradient) * storedDensity_[elementIdx];
unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
unsigned int flowDirectionAxis = asImp_().flowDirectionAxis(elementIdx);
Scalar velocityGradient = asImp_().velocityGradient(elementIdx, flowDirectionAxis, wallNormalAxis);
return mixingLength * mixingLength * abs(velocityGradient) * storedDensity(elementIdx);
}
//! \brief Returns the wall shear stress velocity
......@@ -260,11 +263,11 @@ public:
{
using std::abs;
using std::sqrt;
unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx];
unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementIdx];
unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementIdx];
return sqrt(asImp_().kinematicViscosity_[wallElementIdx]
* abs(asImp_().velocityGradients_[wallElementIdx][flowNormalAxis][wallNormalAxis]));
unsigned int wallElementIdx = asImp_().wallElementIndex(elementIdx);
unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
unsigned int flowDirectionAxis = asImp_().flowDirectionAxis(elementIdx);
return sqrt(asImp_().kinematicViscosity(wallElementIdx)
* abs(asImp_().velocityGradient(wallElementIdx, flowDirectionAxis, wallNormalAxis)));
}
//! \brief Returns the nominal wall shear stress velocity (accounts for poor approximation of viscous sublayer)
......@@ -272,9 +275,8 @@ public:
{
using std::pow;
using std::sqrt;
unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]];
return pow(cMu(), 0.25)
* sqrt(storedTurbulentKineticEnergy_[matchingPointIdx]);
unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx));
return pow(cMu(), 0.25) * sqrt(storedTurbulentKineticEnergy(matchingPointIndex));
}
/*!
......@@ -283,7 +285,7 @@ public:
const Scalar dissipationWallFunction(unsigned int elementIdx) const
{
return uStarNominal(elementIdx) * uStarNominal(elementIdx) * uStarNominal(elementIdx)
/ asImp_().karmanConstant() / asImp_().wallDistance_[elementIdx];
/ asImp_().karmanConstant() / asImp_().wallDistance(elementIdx);
}
/*!
......@@ -291,9 +293,8 @@ public:
*/
const Scalar turbulentKineticEnergyWallFunction(unsigned int elementIdx) const
{
unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx];
unsigned int matchingPointIdx = matchingPointIdx_[wallElementIdx];
return storedTurbulentKineticEnergy_[matchingPointIdx];
unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx));
return storedTurbulentKineticEnergy(matchingPointIndex);
}
//! \brief Returns the nominal wall shear stress (accounts for poor approximation of viscous sublayer)
......@@ -454,21 +455,46 @@ public:
//! \brief Returns the \f$ C_{\mu} \f$ constant
const Scalar cMu() const
{ return 0.09; }
Scalar yPlusThreshold() const
{
return 0.09;
static const Scalar yPlusThreshold = getParamFromGroup<Scalar>(this->paramGroup(), "KEpsilon.YPlusThreshold", 30);
return yPlusThreshold;
}
public:
bool useStoredEddyViscosity() const
{
static const bool useStoredEddyViscosity = getParamFromGroup<bool>(this->paramGroup(), "RANS.UseStoredEddyViscosity", false);
return useStoredEddyViscosity;
}
Scalar storedDensity(const int elementIdx) const
{ return storedDensity_[elementIdx]; }
Scalar storedDissipation(const int elementIdx) const
{ return storedDissipation_[elementIdx]; }
Scalar storedTurbulentKineticEnergy(const int elementIdx) const
{ return storedTurbulentKineticEnergy_[elementIdx]; }
Scalar storedDynamicEddyViscosity(const int elementIdx) const
{ return storedDynamicEddyViscosity_[elementIdx]; }
Scalar zeroEqDynamicEddyViscosity(const int elementIdx) const
{ return zeroEqDynamicEddyViscosity_[elementIdx]; }
unsigned int matchingPointIdx(const int elementIdx) const
{ return matchingPointIdx_[elementIdx]; }
private:
std::vector<unsigned int> matchingPointIdx_;
std::vector<Scalar> storedDensity_;
std::vector<Scalar> storedDissipation_;
std::vector<Scalar> storedTurbulentKineticEnergy_;
std::vector<Scalar> storedDynamicEddyViscosity_;
std::vector<Scalar> zeroEqDynamicEddyViscosity_;
bool useStoredEddyViscosity_;
Scalar yPlusThreshold_;
private:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
......
......@@ -85,21 +85,21 @@ public:
inNearWallRegion_ = problem.inNearWallRegion(RANSParentType::elementIdx());
turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx];
dissipation_ = elemSol[0][Indices::dissipationIdx];
storedDissipation_ = problem.storedDissipation_[RANSParentType::elementIdx()];
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementIdx()];
stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()];
storedDissipation_ = problem.storedDissipation(RANSParentType::elementIdx());
storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy(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 flowDirectionAxis = problem.flowDirectionAxis(RANSParentType::elementIdx());
yPlusNominal_ = RANSParentType::wallDistance() * uStarNominal / problem.kinematicViscosity(RANSParentType::elementIdx());
uPlusNominal_ = RANSParentType::ccVelocityVector()[flowDirectionAxis] / uStarNominal;
cMu_ = problem.cMu();
if (problem.useStoredEddyViscosity_)
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]);
if (problem.useStoredEddyViscosity())
RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity(RANSParentType::elementIdx()));
else
RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity());
if (inNearWallRegion_)
{
RANSParentType::setDynamicEddyViscosity_(problem.zeroEqDynamicEddyViscosity_[RANSParentType::elementIdx()]);
RANSParentType::setDynamicEddyViscosity_(problem.zeroEqDynamicEddyViscosity(RANSParentType::elementIdx()));
}
RANSParentType::calculateEddyDiffusivity(problem);
RANSParentType::calculateEddyThermalConductivity(problem);
......@@ -118,57 +118,43 @@ public:
* \brief Returns the turbulent kinetic energy \f$ m^2/s^2 \f$
*/
Scalar turbulentKineticEnergy() const
{
return turbulentKineticEnergy_;
}
{ return turbulentKineticEnergy_; }
/*!
* \brief Returns an effective dissipation \f$ m^2/s^3 \f$
*/
Scalar dissipation() const
{
return dissipation_;
}
{ return dissipation_; }
/*!
* \brief Returns the turbulent kinetic energy \f$ m^2/s^2 \f$
*/
Scalar storedTurbulentKineticEnergy() const
{
return storedTurbulentKineticEnergy_;
}
{ return storedTurbulentKineticEnergy_; }
/*!
* \brief Returns an effective dissipation \f$ m^2/s^3 \f$
*/
Scalar storedDissipation() const
{
return storedDissipation_;
}
{ return storedDissipation_; }
/*!