Commit f02f564a by Thomas Fetzer

[kepsilon] Laufer's pipe problem works

parent d930b63f
 ... ... @@ -26,55 +26,38 @@ * * The k-epsilon models calculate the eddy viscosity with two additional PDEs, * one for the turbulent kinetic energy (k) and for the dissipation (\f$\varepsilon \f$). * \todo fix everything below here * The model uses the one proposed by Chien \cite Chien1982a. * A good overview and additional models are given in Patel et al. \cite Patel1985a. * * The turbulent kinetic energy balance is identical with the one from the k-epsilon model, * but the dissipation includes a dampening function (\f$D_\varepsilon \f$): * \f$\varepsilon = \tilde{\varepsilon} + D_\varepsilon \f$: * The model uses the one proposed by Launder and Sharma \cite Launder1994a. * * The turbulent kinetic energy balance is: * \f[ * \frac{\partial \left( k \right)}{\partial t} * + \nabla \cdot \left( \textbf{v} k \right) * - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_\text{k}} \right) \nabla k \right) * - 2 \nu_\text{t} \textbf{S} \cdot \textbf{S} * + \tilde{\varepsilon} * + D_\varepsilon * + \varepsilon * = 0 * \f]. * * The dissipation balance is changed by introducing additional functions * (\f$E_\text{k}\f$, \f$f_1 \f$, and \f$f_2 \f$) to account for a dampening towards the wall: * The dissipation balance is: * \f[ * \frac{\partial \left( \tilde{\varepsilon} \right)}{\partial t} * + \nabla \cdot \left( \textbf{v} \tilde{\varepsilon} \right) * - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_{\varepsilon}} \right) \nabla \tilde{\varepsilon} \right) * - C_{1\tilde{\varepsilon}} f_1 \frac{\tilde{\varepsilon}}{k} 2 \nu_\text{t} \textbf{S} \cdot \textbf{S} * + C_{2\tilde{\varepsilon}} f_2 \frac{\tilde{\varepsilon}^2}{k} * - E_\text{k} * \frac{\partial \left( \varepsilon \right)}{\partial t} * + \nabla \cdot \left( \textbf{v} \varepsilon \right) * - \nabla \cdot \left( \left( \nu + \frac{\nu_\text{t}}{\sigma_{\varepsilon}} \right) \nabla \varepsilon \right) * - C_{1\varepsilon} \frac{\varepsilon}{k} 2 \nu_\text{t} \textbf{S} \cdot \textbf{S} * + C_{2\varepsilon} \frac{\varepsilon^2}{k} * = 0 * \f]. * * The kinematic eddy viscosity \f$\nu_\text{t} \f$ is dampened by \f$f_\mu \f$: * The kinematic eddy viscosity \f$\nu_\text{t} \f$ is: * \f[ * \nu_\text{t} = C_\mu f_\mu \frac{k^2}{\tilde{\varepsilon}} * \nu_\text{t} = C_\mu \frac{k^2}{\tilde{\varepsilon}} * \f]. * * The auxiliary and dampening functions are defined as: * \f[ D_\varepsilon = 2 \nu \nicefrac{k}{y^2} \f] * \f[ E_\text{k} = -2 \nu \frac{\tilde{\varepsilon}}{y^2} \exp \left( -0.5 y^+ \right) \f] * \f[ f_1 = 1 \f] * \f[ f_2 = 1 - 0.22 \exp \left( - \left( \frac{\mathit{Re}_\text{t}}{6} \right)^2 \right) \f] * \f[ f_\mu = 1 - \exp \left( -0.0115 y^+ \right) \f] * \f[ \mathit{Re}_\text{t} = \frac{k^2}{\nu \tilde{\varepsilon}} \f] * . * * Finally, the model is closed with the following constants: * \f[ \sigma_\text{k} = 1.00 \f] * \f[ \sigma_\varepsilon =1.30 \f] * \f[ C_{1\tilde{\varepsilon}} = 1.35 \f] * \f[ C_{2\tilde{\varepsilon}} = 1.80 \f] * \f[ C_{1\varepsilon} = 1.44 \f] * \f[ C_{2\varepsilon} = 1.92 \f] * \f[ C_\mu = 0.09 \f] */ ... ...
 ... ... @@ -65,7 +65,8 @@ public: KEpsilonProblem(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) { useStoredEddyViscosity_ = getParamFromGroup(this->paramGroup(), "RANS.UseStoredEddyViscosity", true); yPlusThreshold_ = getParamFromGroup(this->paramGroup(), "KEpsilon.YPlusThreshold", 30); useStoredEddyViscosity_ = getParamFromGroup(this->paramGroup(), "RANS.UseStoredEddyViscosity", false); } /*! ... ... @@ -76,9 +77,12 @@ public: ParentType::updateStaticWallProperties(); // update size and initial values of the global vectors storedDissipationTilde_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); storedDynamicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); matchingPointID_.resize(this->fvGridGeometry().elementMapper().size(), 0); storedDensity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); storedDissipation_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); storedTurbulentKineticEnergy_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); storedDynamicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); zeroEqDynamicEddyViscosity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); } /*! ... ... @@ -90,6 +94,7 @@ public: { ParentType::updateDynamicWallProperties(curSol); // update the stored eddy viscosities for (const auto& element : elements(this->fvGridGeometry().gridView())) { unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); ... ... @@ -103,21 +108,166 @@ public: PrimaryVariables priVars = makePriVarsFromCellCenterPriVars(cellCenterPriVars); auto elemSol = elementSolution(std::move(priVars)); // NOTE: first update the turbulence quantities storedDissipationTilde_[elementID] = elemSol[0][Indices::dissipationEqIdx]; storedDissipation_[elementID] = elemSol[0][Indices::dissipationEqIdx]; storedTurbulentKineticEnergy_[elementID] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; // NOTE: then update the volVars VolumeVariables volVars; volVars.update(elemSol, asImp_(), element, scv); storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity(); storedDensity_[elementID] = volVars.density(); } } // get matching point for k-epsilon wall function unsigned int numElementsInWallRegion = 0; for (const auto& element : elements(this->fvGridGeometry().gridView())) { unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementID]; unsigned int neighborID0 = asImp_().neighborID_[elementID][wallNormalAxis][0]; unsigned int neighborID1 = asImp_().neighborID_[elementID][wallNormalAxis][1]; numElementsInWallRegion = inNearWallRegion(elementID) ? numElementsInWallRegion + 1 : numElementsInWallRegion + 0; if ((inNearWallRegion(elementID) && (inNearWallRegion(neighborID0) != inNearWallRegion(neighborID1))) || (inNearWallRegion(elementID) && (asImp_().wallElementID_[neighborID0] != asImp_().wallElementID_[neighborID1])) || (elementID == asImp_().wallElementID_[elementID] && (!inNearWallRegion(neighborID0) || !inNearWallRegion(neighborID1)))) { matchingPointID_[asImp_().wallElementID_[elementID]] = elementID; } } std::cout << "numElementsInWallRegion: " << numElementsInWallRegion << std::endl; // calculate the potential zeroeq eddy viscosities for two-layer model for (const auto& element : elements(this->fvGridGeometry().gridView())) { unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); zeroEqDynamicEddyViscosity_[elementID] = zeroEqEddyViscosityModel(elementID); } // then make them match at the matching point static const auto enableZeroEqScaling = getParamFromGroup(this->paramGroup(), "KEpsilon.EnableZeroEqScaling", true); if (enableZeroEqScaling) { for (const auto& element : elements(this->fvGridGeometry().gridView())) { unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); unsigned int matchingPointID = matchingPointID_[asImp_().wallElementID_[elementID]]; Scalar scalingFactor = storedDynamicEddyViscosity_[matchingPointID] / zeroEqDynamicEddyViscosity_[matchingPointID]; if (!isMatchingPoint(elementID) && !std::isnan(scalingFactor) && !std::isinf(scalingFactor)) { zeroEqDynamicEddyViscosity_[elementID] *= scalingFactor; } } for (const auto& element : elements(this->fvGridGeometry().gridView())) { unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); unsigned int matchingPointID = matchingPointID_[asImp_().wallElementID_[elementID]]; if (isMatchingPoint(elementID)) { zeroEqDynamicEddyViscosity_[matchingPointID] = storedDynamicEddyViscosity_[matchingPointID]; } } } } /*! * \brief Returns if an element is located in the near-wall region */ const bool inNearWallRegion(unsigned int elementID) const { return yPlus(elementID) < yPlusThreshold_; } /*! * \brief Returns if an element is the matching point */ const bool isMatchingPoint(unsigned int elementID) const { return matchingPointID_[asImp_().wallElementID_[elementID]] == elementID; } /*! * \brief Returns the \f$y^+ \f$ value at an element center */ const Scalar yPlus(unsigned int elementID) const { return asImp_().wallDistance_[elementID] * uStar(elementID) / asImp_().kinematicViscosity_[elementID]; } /*! * \brief Returns the kinematic eddy viscosity of a 0-Eq. model */ const Scalar zeroEqEddyViscosityModel(unsigned int elementID) const { using std::abs; using std::exp; using std::sqrt; // use VanDriest's model Scalar yPlusValue = yPlus(elementID); Scalar mixingLength = 0.0; if (yPlusValue > 0.0) { mixingLength = asImp_().karmanConstant() * asImp_().wallDistance_[elementID] * (1.0 - exp(-yPlusValue / 26.0 )) / sqrt(1.0 - exp(-0.26 * yPlusValue)); } unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementID]; unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementID]; Scalar velocityGradient = asImp_().velocityGradients_[elementID][flowNormalAxis][wallNormalAxis]; return mixingLength * mixingLength * abs(velocityGradient) * storedDensity_[elementID]; } //! \brief Returns the wall shear stress velocity const Scalar uStar(unsigned int elementID) const { using std::abs; using std::sqrt; unsigned int wallElementID = asImp_().wallElementID_[elementID]; unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementID]; unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementID]; return sqrt(asImp_().kinematicViscosity_[wallElementID] * abs(asImp_().velocityGradients_[wallElementID][flowNormalAxis][wallNormalAxis])); } //! \brief Returns the nominal wall shear stress velocity (accounts for poor approximation of viscous sublayer) const Scalar uStarNominal(unsigned int elementID) const { using std::max; using std::pow; using std::sqrt; unsigned int matchingPointID = matchingPointID_[asImp_().wallElementID_[elementID]]; return pow(0.09/*c_mu*/, 0.25) * sqrt(max(storedTurbulentKineticEnergy_[matchingPointID],1e-8)); } /*! * \brief Returns the dissipation calculated from the wall function consideration */ const Scalar dissipationWallFunction(unsigned int elementID) const { return uStarNominal(elementID) * uStarNominal(elementID) * uStarNominal(elementID) / asImp_().karmanConstant() / asImp_().wallDistance_[elementID]; } /*! * \brief Returns the turbulentKineticEnergy calculated from the wall function consideration */ const Scalar turbulentKineticEnergyWallFunction(unsigned int elementID) const { return storedTurbulentKineticEnergy_[elementID]; } public: std::vector storedDissipationTilde_; std::vector storedDynamicEddyViscosity_; std::vector matchingPointID_; std::vector storedDensity_; std::vector storedDissipation_; std::vector storedTurbulentKineticEnergy_; std::vector storedDynamicEddyViscosity_; std::vector zeroEqDynamicEddyViscosity_; bool useStoredEddyViscosity_; Scalar yPlusThreshold_; private: //! Returns the implementation of the problem (i.e. static polymorphism) ... ...
 ... ... @@ -71,10 +71,8 @@ class KEpsilonFluxVariablesImpl(problem.paramGroup(), "KEpsilon.EnableKinematicViscosity_", true); if (kEpsilonEnableKinematicViscosity_) { insideCoeff_k += insideVolVars.kinematicViscosity(); outsideCoeff_k += outsideVolVars.kinematicViscosity(); insideCoeff_e += insideVolVars.kinematicViscosity(); outsideCoeff_e += outsideVolVars.kinematicViscosity(); } // scale by extrusion factor insideCoeff_k *= insideVolVars.extrusionFactor(); ... ... @@ -147,20 +150,24 @@ public: } const auto bcTypes = problem.boundaryTypes(element, scvf); if (!(scvf.boundary() && (bcTypes.isOutflow(turbulentKineticEnergyEqIdx) if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::turbulentKineticEnergyEqIdx) || bcTypes.isSymmetry()))) { flux[turbulentKineticEnergyEqIdx - ModelTraits::dim()] += coeff_k / distance * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) * scvf.area(); if (!insideVolVars.inNearWallRegion() || !insideVolVars.inNearWallRegion()) { flux[turbulentKineticEnergyEqIdx] += coeff_k / distance * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) * scvf.area(); } } if (!(scvf.boundary() && (bcTypes.isOutflow(dissipationEqIdx) if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::dissipationEqIdx) || bcTypes.isSymmetry()))) { flux[dissipationEqIdx - ModelTraits::dim()] flux[dissipationEqIdx] += coeff_e / distance * (insideVolVars.dissipationTilde() - outsideVolVars.dissipationTilde()) * (insideVolVars.dissipation() - outsideVolVars.dissipation()) * scvf.area(); } return flux; ... ...
 ... ... @@ -44,6 +44,7 @@ class KEpsilonResidualImpl; using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); ... ... @@ -69,9 +70,18 @@ class KEpsilonResidualImpl exclude both terms (according to local equilibrium hypothesis, see FLUENT) if (!problem.isMatchingPoint(elementID)) { source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.kinematicEddyViscosity() * volVars.stressTensorScalarProduct(); } source[dissipationEqIdx] += volVars.cOneEpsilon() * dissipation / turbulentKineticEnergy * 2.0 * volVars.kinematicEddyViscosity() * volVars.stressTensorScalarProduct(); // destruction source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde(); source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo() * volVars.dissipationTilde() * volVars.dissipationTilde() / volVars.turbulentKineticEnergy(); // dampening functions source[turbulentKineticEnergyEqIdx] -= volVars.dValue(); source[dissipationEqIdx] += volVars.eValue(); // turbulence production is equal to dissipation -> exclude both terms (according to local equilibrium hypothesis, see FLUENT) if (!problem.isMatchingPoint(elementID)) { source[turbulentKineticEnergyEqIdx] -= dissipation; } source[dissipationEqIdx] -= volVars.cTwoEpsilon() * dissipation * dissipation / turbulentKineticEnergy; return source; } protected: /*! * \brief Evaluate boundary conditions for a cell center dof */ template void evalBoundaryForCellCenter_(CellCenterResidual& residual, const Problem& problem, const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const ElementFaceVariables& elemFaceVars, const ElementBoundaryTypes& elemBcTypes, const ElementFluxVariablesCache& elemFluxVarsCache) const { BaseLocalResidual::evalBoundaryForCellCenter_(residual, problem, element, fvGeometry, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache); for (auto&& scvf : scvfs(fvGeometry)) { unsigned int elementID = problem.fvGridGeometry().elementMapper().index(element); const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); const auto& insideVolVars = elemVolVars[insideScv]; // fixed value for the turbulent kinetic energy if(problem.inNearWallRegion(elementID) && !problem.isMatchingPoint(elementID)) { // unsigned int matchingPointID = problem.matchingPointID_[elementID]; residual[Indices::turbulentKineticEnergyEqIdx - cellCenterOffset] = insideVolVars.turbulentKineticEnergy() - problem.turbulentKineticEnergyWallFunction(elementID); } // fixed value for the dissipation if(problem.inNearWallRegion(elementID)) { residual[Indices::dissipationEqIdx - cellCenterOffset] = insideVolVars.dissipation() - problem.dissipationWallFunction(elementID); } } } }; } ... ...
 ... ... @@ -93,15 +93,21 @@ public: const SubControlVolume& scv) { RANSParentType::updateRANSProperties(elemSol, problem, element, scv); isMatchingPoint_ = problem.isMatchingPoint(RANSParentType::elementID()); inNearWallRegion_ = problem.inNearWallRegion(RANSParentType::elementID()); turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx]; dissipationTilde_ = elemSol[0][Indices::dissipationIdx]; storedDissipationTilde_ = problem.storedDissipationTilde_[RANSParentType::elementID()]; dissipation_ = elemSol[0][Indices::dissipationIdx]; storedDissipation_ = problem.storedDissipation_[RANSParentType::elementID()]; storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()]; stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()]; if (problem.useStoredEddyViscosity_) dynamicEddyViscosity_ = problem.storedDynamicEddyViscosity_[RANSParentType::elementID()]; else dynamicEddyViscosity_ = calculateEddyViscosity(); if (inNearWallRegion_ && !isMatchingPoint_) { dynamicEddyViscosity_ = problem.zeroEqDynamicEddyViscosity_[RANSParentType::elementID()]; } calculateEddyDiffusivity(problem); } ... ... @@ -134,8 +140,8 @@ public: */ Scalar calculateEddyViscosity() { return cMu() * fMu() * turbulentKineticEnergy() * turbulentKineticEnergy() / dissipationTilde() * NavierStokesParentType::density(); return cMu() * turbulentKineticEnergy() * turbulentKineticEnergy() / dissipation() * NavierStokesParentType::density(); } /*! ... ... @@ -162,9 +168,9 @@ public: /*! * \brief Returns an effective dissipation \f$m^2/s^3 \f$ */ Scalar dissipationTilde() const Scalar dissipation() const { return dissipationTilde_; return dissipation_; } /*! ... ... @@ -178,9 +184,9 @@ public: /*! * \brief Returns an effective dissipation \f$m^2/s^3 \f$ */ Scalar storedDissipationTilde() const Scalar storedDissipation() const { return storedDissipationTilde_; return storedDissipation_; } /*! ... ... @@ -191,19 +197,20 @@ public: return stressTensorScalarProduct_; } //! \brief Returns the \$f Re_\textrm{T} \$f value const Scalar reT() const /* * \brief Returns if an element is located in the near-wall region */ bool inNearWallRegion() const { return turbulentKineticEnergy() * turbulentKineticEnergy() / RANSParentType::kinematicViscosity() / dissipationTilde(); return inNearWallRegion_; } //! \brief Returns the \$f Re_\textrm{y} \$f value const Scalar reY() const /*! * \brief Returns if an element is the matching point */ Scalar isMatchingPoint() const { using std::sqrt; return sqrt(turbulentKineticEnergy()) * RANSParentType::wallDistance() / RANSParentType::kinematicViscosity(); return isMatchingPoint_; } //! \brief Returns the \$f C_\mu \$f constant ... ... @@ -218,48 +225,13 @@ public: const Scalar sigmaEpsilon() const { return 1.3; } //! \brief Returns the \$f C_{1\tilde{\varepsilon}} \$f constant //! \brief Returns the \$f C_{1\varepsilon} \$f constant const Scalar