| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup KEpsilonModel | ||
| 10 | * \brief K-epsilon turbulence problem base class. | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_KEPSILON_PROBLEM_HH | ||
| 13 | #define DUMUX_KEPSILON_PROBLEM_HH | ||
| 14 | |||
| 15 | #include <numeric> | ||
| 16 | |||
| 17 | #include <dumux/common/properties.hh> | ||
| 18 | #include <dumux/common/staggeredfvproblem.hh> | ||
| 19 | #include <dumux/discretization/localview.hh> | ||
| 20 | #include <dumux/discretization/staggered/elementsolution.hh> | ||
| 21 | #include <dumux/discretization/method.hh> | ||
| 22 | #include <dumux/freeflow/rans/problem.hh> | ||
| 23 | #include <dumux/freeflow/turbulencemodel.hh> | ||
| 24 | |||
| 25 | #include "model.hh" | ||
| 26 | |||
| 27 | namespace Dumux { | ||
| 28 | |||
| 29 | /*! | ||
| 30 | * \ingroup KEpsilonModel | ||
| 31 | * \brief K-epsilon turbulence problem base class. | ||
| 32 | * | ||
| 33 | * This implements some base functionality for k-epsilon models. | ||
| 34 | */ | ||
| 35 | template<class TypeTag> | ||
| 36 | class RANSProblemImpl<TypeTag, TurbulenceModel::kepsilon> : public RANSProblemBase<TypeTag> | ||
| 37 | { | ||
| 38 | using ParentType = RANSProblemBase<TypeTag>; | ||
| 39 | using Implementation = GetPropType<TypeTag, Properties::Problem>; | ||
| 40 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 41 | |||
| 42 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 43 | using GridView = typename GridGeometry::GridView; | ||
| 44 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 45 | |||
| 46 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 47 | using GridFaceVariables = typename GridVariables::GridFaceVariables; | ||
| 48 | using ElementFaceVariables = typename GridFaceVariables::LocalView; | ||
| 49 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 50 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 51 | |||
| 52 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 53 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 54 | using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; | ||
| 55 | |||
| 56 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 57 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 58 | using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 59 | using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>; | ||
| 60 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 61 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 62 | using Indices = typename ModelTraits::Indices; | ||
| 63 | |||
| 64 | static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance(); | ||
| 65 | static constexpr bool isCompositional = ModelTraits::numFluidComponents() > 1; | ||
| 66 | |||
| 67 | // account for the offset of the cell center privars within the PrimaryVariables container | ||
| 68 | static constexpr auto cellCenterOffset = ModelTraits::numEq() - CellCenterPrimaryVariables::dimension; | ||
| 69 | static_assert(cellCenterOffset == ModelTraits::dim(), "cellCenterOffset must equal dim for staggered NavierStokes"); | ||
| 70 | |||
| 71 | public: | ||
| 72 | |||
| 73 | //! The constructor sets the gravity, if desired by the user. | ||
| 74 | 8 | RANSProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "") | |
| 75 |
2/6✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
48 | : ParentType(gridGeometry, paramGroup) |
| 76 | 8 | { } | |
| 77 | |||
| 78 | /*! | ||
| 79 | * \brief Correct size of the static (solution independent) wall variables | ||
| 80 | */ | ||
| 81 | 8 | void updateStaticWallProperties() | |
| 82 | { | ||
| 83 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
|
8 | if (!ParentType::isFlatWallBounded()) |
| 84 | { | ||
| 85 | ✗ | DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, k-epsilon models should only be used for " | |
| 86 | << " wall bounded flows with flat channel geometries. " | ||
| 87 | << "\n If your geometry is a flat channel, please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); | ||
| 88 | } | ||
| 89 | |||
| 90 | 8 | ParentType::updateStaticWallProperties(); | |
| 91 | // update size and initial values of the global vectors | ||
| 92 | 32 | matchingPointIdx_.resize(this->gridGeometry().elementMapper().size(), 0); | |
| 93 | 32 | storedDissipation_.resize(this->gridGeometry().elementMapper().size(), 0.0); | |
| 94 | 32 | storedTurbulentKineticEnergy_.resize(this->gridGeometry().elementMapper().size(), 0.0); | |
| 95 | 32 | storedDynamicEddyViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0); | |
| 96 | 32 | zeroEqDynamicEddyViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0); | |
| 97 | 8 | } | |
| 98 | |||
| 99 | /*! | ||
| 100 | * \brief Update the dynamic (solution dependent) relations to the walls | ||
| 101 | * | ||
| 102 | * \param curSol The solution vector. | ||
| 103 | */ | ||
| 104 | template<class SolutionVector> | ||
| 105 | 287 | void updateDynamicWallProperties(const SolutionVector& curSol) | |
| 106 | { | ||
| 107 | 287 | ParentType::updateDynamicWallProperties(curSol); | |
| 108 | |||
| 109 | // update the stored eddy viscosities | ||
| 110 |
0/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
574 | auto fvGeometry = localView(this->gridGeometry()); |
| 111 |
1/9✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 64903 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
130093 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 112 | { | ||
| 113 | 193848 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 114 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
64616 | fvGeometry.bindElement(element); |
| 115 |
4/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 64616 times.
✓ Branch 3 taken 64616 times.
✓ Branch 4 taken 64616 times.
✓ Branch 5 taken 64616 times.
|
258464 | for (auto&& scv : scvs(fvGeometry)) |
| 116 | { | ||
| 117 | 64616 | const int dofIdx = scv.dofIndex(); | |
| 118 | 129232 | const auto& cellCenterPriVars = curSol[GridGeometry::cellCenterIdx()][dofIdx]; | |
| 119 | 64616 | PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars); | |
| 120 |
0/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
129232 | auto elemSol = elementSolution<typename GridGeometry::LocalView>(std::move(priVars)); |
| 121 | // NOTE: first update the turbulence quantities | ||
| 122 | 129232 | storedDissipation_[elementIdx] = elemSol[0][Indices::dissipationEqIdx]; | |
| 123 | 129232 | storedTurbulentKineticEnergy_[elementIdx] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; | |
| 124 | // NOTE: then update the volVars | ||
| 125 | 64616 | VolumeVariables volVars; | |
| 126 |
1/2✓ Branch 1 taken 64616 times.
✗ Branch 2 not taken.
|
64616 | volVars.update(elemSol, asImp_(), element, scv); |
| 127 |
3/6✓ Branch 0 taken 64616 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 64616 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 64616 times.
✗ Branch 5 not taken.
|
193848 | storedDynamicEddyViscosity_[elementIdx] = volVars.calculateEddyViscosity(); |
| 128 | } | ||
| 129 | } | ||
| 130 | |||
| 131 | // get matching point for k-epsilon wall function | ||
| 132 | 287 | unsigned int numElementsInNearWallRegion = 0; | |
| 133 |
1/9✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 64903 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
130093 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 134 | { | ||
| 135 | 193848 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 136 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
64616 | unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx); |
| 137 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
64616 | unsigned int neighborIndex0 = asImp_().neighborIndex(elementIdx, wallNormalAxis, 0); |
| 138 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
64616 | unsigned int neighborIndex1 = asImp_().neighborIndex(elementIdx, wallNormalAxis, 1); |
| 139 | 64616 | numElementsInNearWallRegion = inNearWallRegion(elementIdx) | |
| 140 |
2/4✓ Branch 1 taken 14819 times.
✓ Branch 2 taken 49797 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
64616 | ? numElementsInNearWallRegion + 1 |
| 141 | : numElementsInNearWallRegion + 0; | ||
| 142 |
4/11✗ Branch 1 not taken.
✓ Branch 2 taken 47479 times.
✓ Branch 3 taken 2318 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 45786 times.
✓ Branch 6 taken 1693 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
114413 | if ((!inNearWallRegion(elementIdx) && (inNearWallRegion(neighborIndex0) || inNearWallRegion(neighborIndex1))) |
| 143 |
4/9✓ Branch 1 taken 45786 times.
✓ Branch 2 taken 14819 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 43496 times.
✓ Branch 5 taken 2290 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
60605 | || (!inNearWallRegion(elementIdx) && elementIdx == asImp_().wallElementIndex(elementIdx)) |
| 144 |
6/13✓ Branch 0 taken 49797 times.
✓ Branch 1 taken 14819 times.
✓ Branch 3 taken 14819 times.
✓ Branch 4 taken 43496 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 14809 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
122931 | || (inNearWallRegion(elementIdx) && (asImp_().wallElementIndex(neighborIndex0) != asImp_().wallElementIndex(neighborIndex1)))) |
| 145 | { | ||
| 146 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
6311 | matchingPointIdx_[asImp_().wallElementIndex(elementIdx)] = elementIdx; |
| 147 | } | ||
| 148 | } | ||
| 149 |
0/4✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
574 | std::cout << "numElementsInNearWallRegion: " << numElementsInNearWallRegion << std::endl; |
| 150 | |||
| 151 | // calculate the potential zeroeq eddy viscosities for two-layer model | ||
| 152 |
1/9✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 64903 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
130093 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 153 | { | ||
| 154 | 193848 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 155 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
64616 | zeroEqDynamicEddyViscosity_[elementIdx] = zeroEqEddyViscosityModel(elementIdx); |
| 156 | } | ||
| 157 | |||
| 158 | // then make them match at the matching point | ||
| 159 |
3/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 279 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
|
295 | static const auto enableZeroEqScaling |
| 160 |
2/4✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
|
24 | = getParamFromGroup<bool>(this->paramGroup(), "KEpsilon.EnableZeroEqScaling", true); |
| 161 |
1/2✓ Branch 0 taken 287 times.
✗ Branch 1 not taken.
|
287 | if (enableZeroEqScaling) |
| 162 | { | ||
| 163 |
1/9✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 64903 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
130093 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 164 | { | ||
| 165 | 193848 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 166 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
64616 | unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx)); |
| 167 | |||
| 168 | 64616 | Scalar scalingFactor = storedDynamicEddyViscosity(matchingPointIndex) | |
| 169 | 193848 | / zeroEqDynamicEddyViscosity_[matchingPointIndex]; | |
| 170 | 129232 | if (!isMatchingPoint(elementIdx) | |
| 171 |
5/6✓ Branch 0 taken 58294 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 57302 times.
✓ Branch 3 taken 992 times.
✓ Branch 4 taken 57302 times.
✓ Branch 5 taken 992 times.
|
58294 | && !std::isnan(scalingFactor) && !std::isinf(scalingFactor)) |
| 172 | { | ||
| 173 | 114604 | zeroEqDynamicEddyViscosity_[elementIdx] *= scalingFactor; | |
| 174 | } | ||
| 175 | } | ||
| 176 |
1/9✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 64903 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
130380 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 177 | { | ||
| 178 | 193848 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 179 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
64616 | unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx)); |
| 180 | 129232 | if (isMatchingPoint(elementIdx)) | |
| 181 | { | ||
| 182 | 6322 | zeroEqDynamicEddyViscosity_[matchingPointIndex] = storedDynamicEddyViscosity(matchingPointIndex); | |
| 183 | } | ||
| 184 | } | ||
| 185 | } | ||
| 186 | 287 | } | |
| 187 | |||
| 188 | /*! | ||
| 189 | * \brief Returns if an element is located in the near-wall region | ||
| 190 | */ | ||
| 191 | 69242772 | bool inNearWallRegion(unsigned int elementIdx) const | |
| 192 | { | ||
| 193 | 69242772 | unsigned int wallElementIdx = asImp_().wallElementIndex(elementIdx); | |
| 194 | 138485544 | unsigned int matchingPointIndex = matchingPointIdx(wallElementIdx); | |
| 195 | 69242772 | return (wallElementIdx == matchingPointIndex) ? yPlusNominal(elementIdx) < yPlusThreshold() | |
| 196 | 48965108 | : yPlus(elementIdx) < yPlusThreshold(); | |
| 197 | } | ||
| 198 | |||
| 199 | /*! | ||
| 200 | * \brief Returns if an element is the matching point | ||
| 201 | */ | ||
| 202 | bool isMatchingPoint(unsigned int elementIdx) const | ||
| 203 |
0/8✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
59299464 | { return matchingPointIdx(asImp_().wallElementIndex(elementIdx)) == elementIdx; } |
| 204 | |||
| 205 | /*! | ||
| 206 | * \brief Returns the \f$ y^+ \f$ value at an element center | ||
| 207 | */ | ||
| 208 | 49029724 | const Scalar yPlus(unsigned int elementIdx) const | |
| 209 | { | ||
| 210 | 98059448 | return asImp_().wallDistance(elementIdx) * uStar(elementIdx) | |
| 211 | 98059448 | / asImp_().kinematicViscosity(elementIdx); | |
| 212 | } | ||
| 213 | /*! | ||
| 214 | * \brief Returns the nominal \f$ y^+ \f$ value at an element center | ||
| 215 | */ | ||
| 216 | 20978856 | const Scalar yPlusNominal(unsigned int elementIdx) const | |
| 217 | { | ||
| 218 | 41957712 | return asImp_().wallDistance(elementIdx) * uStarNominal(elementIdx) | |
| 219 | 41957712 | / asImp_().kinematicViscosity(elementIdx); | |
| 220 | } | ||
| 221 | |||
| 222 | /*! | ||
| 223 | * \brief Returns the kinematic eddy viscosity of a 0-Eq. model | ||
| 224 | */ | ||
| 225 | 64616 | const Scalar zeroEqEddyViscosityModel(unsigned int elementIdx) const | |
| 226 | { | ||
| 227 | using std::abs; | ||
| 228 | using std::exp; | ||
| 229 | using std::sqrt; | ||
| 230 | |||
| 231 | // use VanDriest's model | ||
| 232 | 64616 | Scalar yPlusValue = yPlus(elementIdx); | |
| 233 | 64616 | Scalar mixingLength = 0.0; | |
| 234 |
1/2✓ Branch 0 taken 64616 times.
✗ Branch 1 not taken.
|
64616 | if (yPlusValue > 0.0) |
| 235 | { | ||
| 236 | 64616 | mixingLength = asImp_().karmanConstant() * asImp_().wallDistance(elementIdx) | |
| 237 | 64616 | * (1.0 - exp(-yPlusValue / 26.0 )) | |
| 238 | 64616 | / sqrt(1.0 - exp(-0.26 * yPlusValue)); | |
| 239 | } | ||
| 240 | |||
| 241 | 64616 | unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx); | |
| 242 | 64616 | unsigned int flowDirectionAxis = asImp_().flowDirectionAxis(elementIdx); | |
| 243 | 129232 | Scalar velocityGradient = asImp_().velocityGradient(elementIdx, flowDirectionAxis, wallNormalAxis); | |
| 244 | 193848 | return mixingLength * mixingLength * abs(velocityGradient) * asImp_().storedDensity(elementIdx); | |
| 245 | } | ||
| 246 | |||
| 247 | //! \brief Returns the wall shear stress velocity | ||
| 248 | 49029724 | const Scalar uStar(unsigned int elementIdx) const | |
| 249 | { | ||
| 250 | using std::abs; | ||
| 251 | using std::sqrt; | ||
| 252 | 49029724 | unsigned int wallElementIdx = asImp_().wallElementIndex(elementIdx); | |
| 253 | 49029724 | unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx); | |
| 254 | 49029724 | unsigned int flowDirectionAxis = asImp_().flowDirectionAxis(elementIdx); | |
| 255 | 98059448 | return sqrt(asImp_().kinematicViscosity(wallElementIdx) | |
| 256 | 245148620 | * abs(asImp_().velocityGradient(wallElementIdx, flowDirectionAxis, wallNormalAxis))); | |
| 257 | } | ||
| 258 | |||
| 259 | //! \brief Returns the nominal wall shear stress velocity (accounts for poor approximation of viscous sublayer) | ||
| 260 | 69351972 | const Scalar uStarNominal(unsigned int elementIdx) const | |
| 261 | { | ||
| 262 | using std::pow; | ||
| 263 | using std::sqrt; | ||
| 264 | 69351972 | unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx)); | |
| 265 | 138703944 | return pow(cMu(), 0.25) * sqrt(storedTurbulentKineticEnergy(matchingPointIndex)); | |
| 266 | } | ||
| 267 | |||
| 268 | /*! | ||
| 269 | * \brief Returns the dissipation calculated from the wall function consideration | ||
| 270 | */ | ||
| 271 | 5632408 | const Scalar dissipationWallFunction(unsigned int elementIdx) const | |
| 272 | { | ||
| 273 | 5632408 | return uStarNominal(elementIdx) * uStarNominal(elementIdx) * uStarNominal(elementIdx) | |
| 274 | 11264816 | / asImp_().karmanConstant() / asImp_().wallDistance(elementIdx); | |
| 275 | } | ||
| 276 | |||
| 277 | /*! | ||
| 278 | * \brief Returns the turbulentKineticEnergy calculated from the wall function consideration | ||
| 279 | */ | ||
| 280 | const Scalar turbulentKineticEnergyWallFunction(unsigned int elementIdx) const | ||
| 281 | { | ||
| 282 | 5632408 | unsigned int matchingPointIndex = matchingPointIdx(asImp_().wallElementIndex(elementIdx)); | |
| 283 | 11264816 | return storedTurbulentKineticEnergy(matchingPointIndex); | |
| 284 | } | ||
| 285 | |||
| 286 | //! \brief Returns the nominal wall shear stress (accounts for poor approximation of viscous sublayer) | ||
| 287 | 508848 | const Scalar tangentialMomentumWallFunction(unsigned int elementIdx, Scalar velocity) const | |
| 288 | { | ||
| 289 | using std::log; | ||
| 290 | 508848 | Scalar velocityNominal = uStarNominal(elementIdx) * (1.0 / asImp_().karmanConstant() * log(yPlusNominal(elementIdx)) + 5.0); | |
| 291 | 508848 | return uStarNominal(elementIdx) * uStarNominal(elementIdx) * velocity / velocityNominal; | |
| 292 | } | ||
| 293 | |||
| 294 | //! \brief Checks whether a wall function should be used | ||
| 295 | 11821545 | bool useWallFunction(const Element& element, | |
| 296 | const SubControlVolumeFace& localSubFace, | ||
| 297 | const int& eqIdx) const | ||
| 298 | { | ||
| 299 | 35464635 | unsigned int elementIdx = asImp_().gridGeometry().elementMapper().index(element); | |
| 300 | 11821545 | auto bcTypes = asImp_().boundaryTypes(element, localSubFace); | |
| 301 | |||
| 302 | return bcTypes.hasWall() | ||
| 303 |
2/2✓ Branch 0 taken 1033654 times.
✓ Branch 1 taken 3458027 times.
|
4491681 | && bcTypes.isDirichlet(eqIdx) |
| 304 |
2/2✓ Branch 0 taken 4491681 times.
✓ Branch 1 taken 7329864 times.
|
15279572 | && isMatchingPoint(elementIdx); |
| 305 | } | ||
| 306 | |||
| 307 | //! \brief Returns an additional wall function momentum flux (only needed for RANS models) | ||
| 308 | ✗ | FacePrimaryVariables wallFunction(const Element& element, | |
| 309 | const FVElementGeometry& fvGeometry, | ||
| 310 | const ElementVolumeVariables& elemVolVars, | ||
| 311 | const ElementFaceVariables& elemFaceVars, | ||
| 312 | const SubControlVolumeFace& scvf, | ||
| 313 | const SubControlVolumeFace& localSubFace) const | ||
| 314 | { | ||
| 315 | ✗ | unsigned int elementIdx = asImp_().gridGeometry().elementMapper().index(element); | |
| 316 | ✗ | return FacePrimaryVariables(asImp_().tangentialMomentumWallFunction(elementIdx, elemFaceVars[scvf].velocitySelf()) | |
| 317 | ✗ | * asImp_().storedDensity(elementIdx) ); | |
| 318 | } | ||
| 319 | |||
| 320 | //! \brief Returns the flux for non-isothermal and compositional RANS models | ||
| 321 | template<bool eB = enableEnergyBalance, bool compositional = isCompositional, | ||
| 322 | typename std::enable_if_t<eB && compositional, int> = 0> | ||
| 323 | 73320 | CellCenterPrimaryVariables wallFunction(const Element& element, | |
| 324 | const FVElementGeometry& fvGeometry, | ||
| 325 | const ElementVolumeVariables& elemVolVars, | ||
| 326 | const ElementFaceVariables& elemFaceVars, | ||
| 327 | const SubControlVolumeFace& scvf) const | ||
| 328 | { | ||
| 329 | return wallFunctionComponent(element, fvGeometry, elemVolVars, elemFaceVars, scvf) | ||
| 330 | 73320 | + wallFunctionEnergy(element, fvGeometry, elemVolVars, elemFaceVars, scvf); | |
| 331 | } | ||
| 332 | |||
| 333 | //! \brief Returns the flux for isothermal and compositional RANS models | ||
| 334 | template<bool eB = enableEnergyBalance, bool compositional = isCompositional, | ||
| 335 | typename std::enable_if_t<!eB && compositional, int> = 0> | ||
| 336 | CellCenterPrimaryVariables wallFunction(const Element& element, | ||
| 337 | const FVElementGeometry& fvGeometry, | ||
| 338 | const ElementVolumeVariables& elemVolVars, | ||
| 339 | const ElementFaceVariables& elemFaceVars, | ||
| 340 | const SubControlVolumeFace& scvf) const | ||
| 341 | 45704 | { return wallFunctionComponent(element, fvGeometry, elemVolVars, elemFaceVars, scvf); } | |
| 342 | |||
| 343 | //! \brief Returns the flux for non-isothermal RANS models | ||
| 344 | template<bool eB = enableEnergyBalance, bool compositional = isCompositional, | ||
| 345 | typename std::enable_if_t<eB && !compositional, int> = 0> | ||
| 346 | CellCenterPrimaryVariables wallFunction(const Element& element, | ||
| 347 | const FVElementGeometry& fvGeometry, | ||
| 348 | const ElementVolumeVariables& elemVolVars, | ||
| 349 | const ElementFaceVariables& elemFaceVars, | ||
| 350 | const SubControlVolumeFace& scvf) const | ||
| 351 | ✗ | { return wallFunctionEnergy(element, fvGeometry, elemVolVars, elemFaceVars, scvf); } | |
| 352 | |||
| 353 | //! \brief Returns the flux for isothermal RANS models | ||
| 354 | template<bool eB = enableEnergyBalance, bool compositional = isCompositional, | ||
| 355 | typename std::enable_if_t<!eB && !compositional, int> = 0> | ||
| 356 | ✗ | CellCenterPrimaryVariables wallFunction(const Element& element, | |
| 357 | const FVElementGeometry& fvGeometry, | ||
| 358 | const ElementVolumeVariables& elemVolVars, | ||
| 359 | const ElementFaceVariables& elemFaceVars, | ||
| 360 | const SubControlVolumeFace& scvf) const | ||
| 361 | 1086168 | { return CellCenterPrimaryVariables(0.0); } | |
| 362 | |||
| 363 | //! \brief Returns the component wall-function flux | ||
| 364 | 119024 | CellCenterPrimaryVariables wallFunctionComponent(const Element& element, | |
| 365 | const FVElementGeometry& fvGeometry, | ||
| 366 | const ElementVolumeVariables& elemVolVars, | ||
| 367 | const ElementFaceVariables& elemFaceVars, | ||
| 368 | const SubControlVolumeFace& scvf) const | ||
| 369 | { | ||
| 370 | using std::log; | ||
| 371 | 119024 | auto wallFunctionFlux = CellCenterPrimaryVariables(0.0); | |
| 372 | 238048 | unsigned int elementIdx = asImp_().gridGeometry().elementMapper().index(element); | |
| 373 | |||
| 374 | // component mass fluxes | ||
| 375 |
2/2✓ Branch 1 taken 238048 times.
✓ Branch 2 taken 119024 times.
|
357072 | for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++compIdx) |
| 376 | { | ||
| 377 |
2/2✓ Branch 0 taken 119024 times.
✓ Branch 1 taken 119024 times.
|
238048 | if (ModelTraits::replaceCompEqIdx() == compIdx) |
| 378 | continue; | ||
| 379 | |||
| 380 | 238048 | Scalar schmidtNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity() | |
| 381 | 357072 | / elemVolVars[scvf.insideScvIdx()].diffusionCoefficient(0, 0, compIdx); | |
| 382 | 119024 | Scalar moleToMassConversionFactor = ModelTraits::useMoles() | |
| 383 | ? 1.0 : FluidSystem::molarMass(compIdx); | ||
| 384 | 238048 | wallFunctionFlux[compIdx] += | |
| 385 | 238048 | -1.0 * (asImp_().dirichlet(element, scvf)[Indices::conti0EqIdx + compIdx] | |
| 386 | 238048 | - elemVolVars[scvf.insideScvIdx()].moleFraction(compIdx)) | |
| 387 | 238048 | * elemVolVars[scvf.insideScvIdx()].molarDensity() | |
| 388 | * moleToMassConversionFactor | ||
| 389 | 119024 | * uStarNominal(elementIdx) | |
| 390 | 119024 | / asImp_().turbulentSchmidtNumber() | |
| 391 | 119024 | / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementIdx) * 9.793) | |
| 392 | 119024 | + pFunction(schmidtNumber, asImp_().turbulentSchmidtNumber())); | |
| 393 | } | ||
| 394 | |||
| 395 | if (ModelTraits::replaceCompEqIdx() < ModelTraits::numFluidComponents()) | ||
| 396 | { | ||
| 397 | 238048 | wallFunctionFlux[ModelTraits::replaceCompEqIdx()] = | |
| 398 | 119024 | -std::accumulate(wallFunctionFlux.begin(), wallFunctionFlux.end(), 0.0); | |
| 399 | } | ||
| 400 | |||
| 401 | 119024 | return wallFunctionFlux; | |
| 402 | } | ||
| 403 | |||
| 404 | //! \brief Returns the energy wall-function flux | ||
| 405 | 73320 | CellCenterPrimaryVariables wallFunctionEnergy(const Element& element, | |
| 406 | const FVElementGeometry& fvGeometry, | ||
| 407 | const ElementVolumeVariables& elemVolVars, | ||
| 408 | const ElementFaceVariables& elemFaceVars, | ||
| 409 | const SubControlVolumeFace& scvf) const | ||
| 410 | { | ||
| 411 | using std::log; | ||
| 412 | 73320 | auto wallFunctionFlux = CellCenterPrimaryVariables(0.0); | |
| 413 | 219960 | unsigned int elementIdx = asImp_().gridGeometry().elementMapper().index(element); | |
| 414 | // energy fluxes | ||
| 415 | 146640 | Scalar prandtlNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity() | |
| 416 | 219960 | * elemVolVars[scvf.insideScvIdx()].density() | |
| 417 | 146640 | * elemVolVars[scvf.insideScvIdx()].heatCapacity() | |
| 418 | 146640 | / elemVolVars[scvf.insideScvIdx()].thermalConductivity(); | |
| 419 | 146640 | wallFunctionFlux[Indices::energyEqIdx - cellCenterOffset] += | |
| 420 | 73320 | -1.0 * (asImp_().dirichlet(element, scvf)[Indices::temperatureIdx] | |
| 421 | 146640 | - elemVolVars[scvf.insideScvIdx()].temperature()) | |
| 422 | 146640 | * elemVolVars[scvf.insideScvIdx()].density() | |
| 423 | 146640 | * elemVolVars[scvf.insideScvIdx()].heatCapacity() | |
| 424 | 73320 | * uStarNominal(elementIdx) | |
| 425 | 73320 | / asImp_().turbulentPrandtlNumber() | |
| 426 | 73320 | / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementIdx) * 9.793) | |
| 427 | 73320 | + pFunction(prandtlNumber, asImp_().turbulentPrandtlNumber())); | |
| 428 | |||
| 429 | 73320 | return wallFunctionFlux; | |
| 430 | } | ||
| 431 | |||
| 432 | //! \brief Returns the value of the P-function after Jayatilleke \cite Versteeg2009a | ||
| 433 | ✗ | const Scalar pFunction(Scalar molecularNumber, Scalar turbulentNumber) const | |
| 434 | { | ||
| 435 | using std::pow; | ||
| 436 | using std::exp; | ||
| 437 | return 9.24 | ||
| 438 | ✗ | * (pow(molecularNumber / turbulentNumber, 0.75) - 1.0) | |
| 439 | ✗ | * (1.0 + 0.28 * exp(-0.007 * molecularNumber / turbulentNumber)); | |
| 440 | } | ||
| 441 | |||
| 442 | //! \brief Returns the \f$ C_{\mu} \f$ constant | ||
| 443 | ✗ | const Scalar cMu() const | |
| 444 | ✗ | { return 0.09; } | |
| 445 | |||
| 446 | 69242772 | Scalar yPlusThreshold() const | |
| 447 | { | ||
| 448 |
5/8✓ Branch 0 taken 8 times.
✓ Branch 1 taken 69242764 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
|
69242772 | static const Scalar yPlusThreshold = getParamFromGroup<Scalar>(this->paramGroup(), "KEpsilon.YPlusThreshold", 30.0); |
| 449 | 69242772 | return yPlusThreshold; | |
| 450 | } | ||
| 451 | |||
| 452 | 29757004 | bool useStoredEddyViscosity() const | |
| 453 | { | ||
| 454 |
5/8✓ Branch 0 taken 8 times.
✓ Branch 1 taken 29756996 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
|
29757004 | static const bool useStoredEddyViscosity = getParamFromGroup<bool>(this->paramGroup(), "RANS.UseStoredEddyViscosity", false); |
| 455 | 29757004 | return useStoredEddyViscosity; | |
| 456 | } | ||
| 457 | |||
| 458 | Scalar storedDissipation(const int elementIdx) const | ||
| 459 | 59514008 | { return storedDissipation_[elementIdx]; } | |
| 460 | |||
| 461 | Scalar storedTurbulentKineticEnergy(const int elementIdx) const | ||
| 462 | 65146416 | { return storedTurbulentKineticEnergy_[elementIdx]; } | |
| 463 | |||
| 464 | Scalar storedDynamicEddyViscosity(const int elementIdx) const | ||
| 465 |
0/8✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
77260 | { return storedDynamicEddyViscosity_[elementIdx]; } |
| 466 | |||
| 467 | Scalar zeroEqDynamicEddyViscosity(const int elementIdx) const | ||
| 468 | 13032410 | { return zeroEqDynamicEddyViscosity_[elementIdx]; } | |
| 469 | |||
| 470 | unsigned int matchingPointIdx(const int elementIdx) const | ||
| 471 |
20/88✓ Branch 0 taken 2287071 times.
✓ Branch 1 taken 1170956 times.
✓ Branch 2 taken 2287071 times.
✓ Branch 3 taken 1170956 times.
✓ Branch 6 taken 3543418 times.
✓ Branch 7 taken 25934426 times.
✓ Branch 8 taken 3543418 times.
✓ Branch 9 taken 25934426 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 20277664 times.
✓ Branch 13 taken 48965108 times.
✓ Branch 14 taken 20277664 times.
✓ Branch 15 taken 48965108 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 58294 times.
✓ Branch 19 taken 6322 times.
✓ Branch 20 taken 58294 times.
✓ Branch 21 taken 6322 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 6322 times.
✓ Branch 25 taken 58294 times.
✓ Branch 26 taken 6322 times.
✓ Branch 27 taken 58294 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
|
275394574 | { return matchingPointIdx_[elementIdx]; } |
| 472 | |||
| 473 | private: | ||
| 474 | std::vector<unsigned int> matchingPointIdx_; | ||
| 475 | std::vector<Scalar> storedDissipation_; | ||
| 476 | std::vector<Scalar> storedTurbulentKineticEnergy_; | ||
| 477 | std::vector<Scalar> storedDynamicEddyViscosity_; | ||
| 478 | std::vector<Scalar> zeroEqDynamicEddyViscosity_; | ||
| 479 | |||
| 480 | //! Returns the implementation of the problem (i.e. static polymorphism) | ||
| 481 | Implementation &asImp_() | ||
| 482 | { return *static_cast<Implementation *>(this); } | ||
| 483 | |||
| 484 | //! \copydoc asImp_() | ||
| 485 | const Implementation &asImp_() const | ||
| 486 | { return *static_cast<const Implementation *>(this); } | ||
| 487 | }; | ||
| 488 | |||
| 489 | } // end namespace Dumux | ||
| 490 | |||
| 491 | #endif | ||
| 492 |