| 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 ZeroEqModel | ||
| 10 | * \brief Zero-equation turbulence problem base class. | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_ZEROEQ_PROBLEM_HH | ||
| 13 | #define DUMUX_ZEROEQ_PROBLEM_HH | ||
| 14 | |||
| 15 | #include <string> | ||
| 16 | |||
| 17 | #include <dune/common/math.hh> | ||
| 18 | #include <dune/common/stdstreams.hh> | ||
| 19 | |||
| 20 | #include <dumux/common/properties.hh> | ||
| 21 | #include <dumux/common/staggeredfvproblem.hh> | ||
| 22 | #include <dumux/discretization/localview.hh> | ||
| 23 | #include <dumux/discretization/staggered/elementsolution.hh> | ||
| 24 | #include <dumux/discretization/method.hh> | ||
| 25 | #include <dumux/freeflow/rans/problem.hh> | ||
| 26 | #include <dumux/freeflow/turbulencemodel.hh> | ||
| 27 | |||
| 28 | #include "model.hh" | ||
| 29 | |||
| 30 | namespace Dumux { | ||
| 31 | |||
| 32 | /*! | ||
| 33 | * \ingroup ZeroEqModel | ||
| 34 | * \brief Zero-equation turbulence problem base class. | ||
| 35 | * | ||
| 36 | * This implements some base functionality for zero-equation models | ||
| 37 | * and a routine for the determining the eddy viscosity of the Baldwin-Lomax model. | ||
| 38 | */ | ||
| 39 | template<class TypeTag> | ||
| 40 | class RANSProblemImpl<TypeTag, TurbulenceModel::zeroeq> : public RANSProblemBase<TypeTag> | ||
| 41 | { | ||
| 42 | using ParentType = RANSProblemBase<TypeTag>; | ||
| 43 | using Implementation = GetPropType<TypeTag, Properties::Problem>; | ||
| 44 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 45 | using Grid = typename GridView::Grid; | ||
| 46 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 47 | |||
| 48 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 49 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 50 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 51 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 52 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 53 | using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 54 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 55 | |||
| 56 | enum { | ||
| 57 | dim = Grid::dimension, | ||
| 58 | }; | ||
| 59 | using DimVector = Dune::FieldVector<Scalar, dim>; | ||
| 60 | using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>; | ||
| 61 | |||
| 62 | public: | ||
| 63 | /*! | ||
| 64 | * \brief The constructor | ||
| 65 | * \param gridGeometry The finite volume grid geometry | ||
| 66 | * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") | ||
| 67 | */ | ||
| 68 | 12 | RANSProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "") | |
| 69 |
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.
|
36 | : ParentType(gridGeometry, paramGroup) |
| 70 | 12 | {} | |
| 71 | |||
| 72 | /*! | ||
| 73 | * \brief Correct size of the static (solution independent) wall variables | ||
| 74 | */ | ||
| 75 | 12 | void updateStaticWallProperties() | |
| 76 | { | ||
| 77 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
|
12 | if (!ParentType::isFlatWallBounded()) |
| 78 | { | ||
| 79 | ✗ | DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, zero-eq models should only be used for flat channel geometries. " | |
| 80 | << "\n If your geometry is a flat channel, please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); | ||
| 81 | } | ||
| 82 | |||
| 83 | 12 | ParentType::updateStaticWallProperties(); | |
| 84 | |||
| 85 | // update size and initial values of the global vectors | ||
| 86 | 48 | kinematicEddyViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0); | |
| 87 | 48 | additionalRoughnessLength_.resize(this->gridGeometry().elementMapper().size(), 0.0); | |
| 88 | 12 | } | |
| 89 | |||
| 90 | /*! | ||
| 91 | * \brief Update the dynamic (solution dependent) relations to the walls | ||
| 92 | * | ||
| 93 | * This calculates the roughness related properties | ||
| 94 | * | ||
| 95 | * \param curSol The solution vector. | ||
| 96 | */ | ||
| 97 | template<class SolutionVector> | ||
| 98 | 139 | void updateDynamicWallProperties(const SolutionVector& curSol) | |
| 99 | { | ||
| 100 | 139 | ParentType::updateDynamicWallProperties(curSol); | |
| 101 | |||
| 102 | // correct roughness lengths if a sand grain roughness is specified | ||
| 103 |
6/14✓ Branch 1 taken 135 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 135 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 135 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 135 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 135 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 135 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
695 | if (hasParam("RANS.SandGrainRoughness")) |
| 104 | ✗ | calculateRoughnessLength_(curSol); | |
| 105 | |||
| 106 | // update routine for specific models | ||
| 107 |
3/4✗ Branch 2 not taken.
✓ Branch 3 taken 135 times.
✓ Branch 4 taken 70 times.
✓ Branch 5 taken 65 times.
|
139 | if (eddyViscosityModel().compare("baldwinLomax") == 0) |
| 108 | 70 | updateBaldwinLomaxProperties(); | |
| 109 | 139 | } | |
| 110 | |||
| 111 | /*! | ||
| 112 | * \brief Update the relations and coefficients for the Baldwin-Lomax turbulence model | ||
| 113 | */ | ||
| 114 | 70 | void updateBaldwinLomaxProperties() | |
| 115 | { | ||
| 116 |
1/2✓ Branch 6 taken 70 times.
✗ Branch 7 not taken.
|
350 | std::vector<Scalar> kinematicEddyViscosityInner(this->gridGeometry().elementMapper().size(), 0.0); |
| 117 |
7/16✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 70 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 70 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 70 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 70 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
420 | std::vector<Scalar> kinematicEddyViscosityOuter(this->gridGeometry().elementMapper().size(), 0.0); |
| 118 |
7/16✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 70 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 70 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 70 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 70 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
420 | std::vector<Scalar> kinematicEddyViscosityDifference(this->gridGeometry().elementMapper().size(), 0.0); |
| 119 |
7/18✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 70 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 70 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 70 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 70 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
420 | std::vector<Scalar> switchingPosition(this->gridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max()); |
| 120 | |||
| 121 | using std::abs; | ||
| 122 | using std::exp; | ||
| 123 | using std::min; | ||
| 124 | using std::sqrt; | ||
| 125 | using Dune::power; | ||
| 126 | 70 | const Scalar aPlus = 26.0; | |
| 127 | 70 | const Scalar k = 0.0168; | |
| 128 | 70 | const Scalar cCP = 1.6; | |
| 129 | 70 | const Scalar cWake = 0.25; | |
| 130 | 70 | const Scalar cKleb = 0.3; | |
| 131 | |||
| 132 |
2/6✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
140 | std::vector<Scalar> storedFMax; |
| 133 |
2/4✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70 times.
✗ Branch 4 not taken.
|
140 | std::vector<Scalar> storedYFMax; |
| 134 |
4/8✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 70 times.
✗ Branch 11 not taken.
|
280 | storedFMax.resize(this->gridGeometry().elementMapper().size(), 0.0); |
| 135 |
4/8✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 70 times.
✗ Branch 11 not taken.
|
280 | storedYFMax.resize(this->gridGeometry().elementMapper().size(), 0.0); |
| 136 | |||
| 137 | // (1) calculate inner viscosity and Klebanoff function | ||
| 138 |
4/8✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 16902 times.
✗ Branch 10 not taken.
|
33874 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 139 | { | ||
| 140 | 50496 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 141 |
1/2✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
|
16832 | Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx); |
| 142 | 33664 | unsigned int flowDirectionAxis = this->flowDirectionAxis(elementIdx); | |
| 143 |
1/2✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
|
16832 | unsigned int wallNormalAxis = this->wallNormalAxis(elementIdx); |
| 144 | |||
| 145 | 33664 | Scalar omegaAbs = abs(this->velocityGradient(elementIdx, flowDirectionAxis, wallNormalAxis) | |
| 146 |
3/6✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16832 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16832 times.
✗ Branch 8 not taken.
|
50496 | - this->velocityGradient(elementIdx, wallNormalAxis, flowDirectionAxis)); |
| 147 | 33664 | Scalar uStar = sqrt(this->kinematicViscosity(asImp_().wallElementIndex(elementIdx)) | |
| 148 |
5/10✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16832 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16832 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 16832 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 16832 times.
✗ Branch 14 not taken.
|
16832 | * abs(this->velocityGradient(asImp_().wallElementIndex(elementIdx), flowDirectionAxis, wallNormalAxis))); |
| 149 |
1/2✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
|
16832 | Scalar yPlus = effectiveWallDistance * uStar / this->kinematicViscosity(elementIdx); |
| 150 | 16832 | Scalar mixingLength = this->karmanConstant() * effectiveWallDistance * (1.0 - exp(-yPlus / aPlus)); | |
| 151 |
1/2✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
|
16832 | kinematicEddyViscosityInner[elementIdx] = mixingLength * mixingLength * omegaAbs; |
| 152 | |||
| 153 | 16832 | Scalar f = effectiveWallDistance * omegaAbs * (1.0 - exp(-yPlus / aPlus)); | |
| 154 |
5/6✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8408 times.
✓ Branch 4 taken 8424 times.
✓ Branch 5 taken 8408 times.
✓ Branch 6 taken 8424 times.
|
16832 | if (f > storedFMax[asImp_().wallElementIndex(elementIdx)]) |
| 155 | { | ||
| 156 |
2/4✓ Branch 1 taken 8408 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8408 times.
✗ Branch 5 not taken.
|
8408 | storedFMax[asImp_().wallElementIndex(elementIdx)] = f; |
| 157 |
1/2✓ Branch 1 taken 8408 times.
✗ Branch 2 not taken.
|
8408 | storedYFMax[asImp_().wallElementIndex(elementIdx)] = effectiveWallDistance; |
| 158 | } | ||
| 159 | } | ||
| 160 | |||
| 161 | // (2) calculate outer viscosity | ||
| 162 |
4/8✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 16902 times.
✗ Branch 10 not taken.
|
33874 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 163 | { | ||
| 164 | 50496 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 165 | 16832 | Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx); | |
| 166 | |||
| 167 | 16832 | Scalar maxVelocityNorm = 0.0; | |
| 168 | 16832 | Scalar minVelocityNorm = 0.0; | |
| 169 |
2/2✓ Branch 0 taken 33664 times.
✓ Branch 1 taken 16832 times.
|
50496 | for (unsigned axisIdx = 0; axisIdx < dim; ++axisIdx) |
| 170 | { | ||
| 171 |
3/6✓ Branch 1 taken 33664 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33664 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 33664 times.
✗ Branch 8 not taken.
|
33664 | maxVelocityNorm += asImp_().velocityMaximum(asImp_().wallElementIndex(elementIdx))[axisIdx] |
| 172 |
3/6✓ Branch 1 taken 33664 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33664 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 33664 times.
✗ Branch 8 not taken.
|
33664 | * asImp_().velocityMaximum(asImp_().wallElementIndex(elementIdx))[axisIdx]; |
| 173 |
3/8✓ Branch 1 taken 33664 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33664 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 33664 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
33664 | minVelocityNorm += asImp_().velocityMinimum(asImp_().wallElementIndex(elementIdx))[axisIdx] |
| 174 |
1/2✓ Branch 1 taken 33664 times.
✗ Branch 2 not taken.
|
33664 | * asImp_().velocityMinimum(asImp_().wallElementIndex(elementIdx))[axisIdx]; |
| 175 | } | ||
| 176 | |||
| 177 | 16832 | Scalar deltaU = sqrt(maxVelocityNorm) - sqrt(minVelocityNorm); | |
| 178 |
2/4✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16832 times.
✗ Branch 5 not taken.
|
16832 | Scalar yFMax = storedYFMax[asImp_().wallElementIndex(elementIdx)]; |
| 179 |
3/4✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2592 times.
✓ Branch 4 taken 14240 times.
|
16832 | Scalar fMax = storedFMax[asImp_().wallElementIndex(elementIdx)]; |
| 180 |
2/2✓ Branch 0 taken 2592 times.
✓ Branch 1 taken 14240 times.
|
16832 | Scalar fWake = min(yFMax * fMax, cWake * yFMax * deltaU * deltaU / fMax); |
| 181 | 16832 | Scalar fKleb = 1.0 / (1.0 + 5.5 * power(cKleb * effectiveWallDistance / yFMax, 6)); | |
| 182 | 16832 | kinematicEddyViscosityOuter[elementIdx] = k * cCP * fWake * fKleb; | |
| 183 | |||
| 184 | 16832 | kinematicEddyViscosityDifference[elementIdx] | |
| 185 | 33664 | = kinematicEddyViscosityInner[elementIdx] - kinematicEddyViscosityOuter[elementIdx]; | |
| 186 | } | ||
| 187 | |||
| 188 | // (3) switching point | ||
| 189 |
4/8✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 16902 times.
✗ Branch 10 not taken.
|
33874 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 190 | { | ||
| 191 | 50496 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 192 |
1/2✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
|
16832 | Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx); |
| 193 | |||
| 194 | // checks if sign switches, by multiplication | ||
| 195 |
4/4✓ Branch 0 taken 9304 times.
✓ Branch 1 taken 7528 times.
✓ Branch 2 taken 9304 times.
✓ Branch 3 taken 7528 times.
|
33664 | Scalar check = kinematicEddyViscosityDifference[asImp_().wallElementIndex(elementIdx)] * kinematicEddyViscosityDifference[elementIdx]; |
| 196 | if (check < 0 // means sign has switched | ||
| 197 |
7/8✓ Branch 0 taken 9304 times.
✓ Branch 1 taken 7528 times.
✓ Branch 3 taken 9304 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 5398 times.
✓ Branch 6 taken 3906 times.
✓ Branch 7 taken 5398 times.
✓ Branch 8 taken 3906 times.
|
16832 | && switchingPosition[asImp_().wallElementIndex(elementIdx)] > effectiveWallDistance) |
| 198 | { | ||
| 199 |
1/2✓ Branch 1 taken 5398 times.
✗ Branch 2 not taken.
|
5398 | switchingPosition[asImp_().wallElementIndex(elementIdx)] = effectiveWallDistance; |
| 200 | } | ||
| 201 | } | ||
| 202 | |||
| 203 | // (4) finally determine eddy viscosity | ||
| 204 |
5/10✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 16902 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 70 times.
✗ Branch 13 not taken.
|
33944 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 205 | { | ||
| 206 | 50496 | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 207 |
1/2✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
|
16832 | Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx); |
| 208 | |||
| 209 |
1/2✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
|
16832 | kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityInner[elementIdx]; |
| 210 |
5/6✓ Branch 1 taken 16832 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10356 times.
✓ Branch 4 taken 6476 times.
✓ Branch 5 taken 10356 times.
✓ Branch 6 taken 6476 times.
|
16832 | if (effectiveWallDistance >= switchingPosition[asImp_().wallElementIndex(elementIdx)]) |
| 211 | { | ||
| 212 | 31068 | kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityOuter[elementIdx]; | |
| 213 | } | ||
| 214 | } | ||
| 215 | 70 | } | |
| 216 | |||
| 217 | 5262887 | std::string eddyViscosityModel() const | |
| 218 | { | ||
| 219 |
5/8✓ Branch 0 taken 8 times.
✓ Branch 1 taken 5262275 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.
|
5262887 | static const std::string eddyViscosityModel = getParamFromGroup<std::string>(this->paramGroup(), "RANS.EddyViscosityModel", "vanDriest"); |
| 220 | 5262887 | return eddyViscosityModel; | |
| 221 | } | ||
| 222 | |||
| 223 | int additionalRoughnessLength(const int elementIdx) const | ||
| 224 |
6/52✓ Branch 3 taken 16832 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 16832 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 16832 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 16832 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 16832 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 16832 times.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 39 not taken.
✗ Branch 40 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 51 not taken.
✗ Branch 52 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 69 not taken.
✗ Branch 70 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.
|
10658952 | { return additionalRoughnessLength_[elementIdx]; } |
| 225 | |||
| 226 | Scalar kinematicEddyViscosity(const int elementIdx) const | ||
| 227 | 5290256 | { return kinematicEddyViscosity_[elementIdx]; } | |
| 228 | |||
| 229 | private: | ||
| 230 | |||
| 231 | template<class SolutionVector> | ||
| 232 | ✗ | void calculateRoughnessLength_(const SolutionVector& curSol) | |
| 233 | { | ||
| 234 | ✗ | bool printedRangeWarning = false; | |
| 235 | ✗ | auto fvGeometry = localView(this->gridGeometry()); | |
| 236 | ✗ | for (const auto& element : elements(this->gridGeometry().gridView())) | |
| 237 | { | ||
| 238 | ✗ | static const Scalar sandGrainRoughness = getParamFromGroup<Scalar>(this->paramGroup(), "RANS.SandGrainRoughness"); | |
| 239 | ✗ | unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); | |
| 240 | ✗ | fvGeometry.bindElement(element); | |
| 241 | |||
| 242 | ✗ | for (auto&& scv : scvs(fvGeometry)) | |
| 243 | { | ||
| 244 | using std::sqrt; | ||
| 245 | using std::exp; | ||
| 246 | |||
| 247 | ✗ | const int dofIdx = scv.dofIndex(); | |
| 248 | |||
| 249 | // construct a privars object from the cell center solution vector | ||
| 250 | ✗ | const auto& cellCenterPriVars = curSol[GridGeometry::cellCenterIdx()][dofIdx]; | |
| 251 | ✗ | PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars); | |
| 252 | ✗ | auto elemSol = elementSolution<typename GridGeometry::LocalView>(std::move(priVars)); | |
| 253 | |||
| 254 | ✗ | VolumeVariables volVars; | |
| 255 | ✗ | volVars.update(elemSol, asImp_(), element, scv); | |
| 256 | |||
| 257 | ✗ | Scalar ksPlus = sandGrainRoughness * volVars.uStar() / volVars.kinematicViscosity(); | |
| 258 | ✗ | if (ksPlus > 0 && eddyViscosityModel().compare("baldwinLomax") == 0) | |
| 259 | { | ||
| 260 | ✗ | DUNE_THROW(Dune::NotImplemented, "Roughness is not implemented for the Baldwin-Lomax model."); | |
| 261 | } | ||
| 262 | ✗ | if (ksPlus > 2000.) | |
| 263 | { | ||
| 264 | ✗ | std::cout << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx)) | |
| 265 | << " is not in the valid range (ksPlus < 2000)," | ||
| 266 | ✗ | << " for high ksPlus values the roughness function reaches a turning point."<< std::endl; | |
| 267 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Unphysical roughness behavior."); | |
| 268 | } | ||
| 269 | ✗ | else if (ksPlus > 0.0 && ksPlus < 4.535 && !printedRangeWarning) | |
| 270 | { | ||
| 271 | ✗ | Dune::dinfo << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx)) | |
| 272 | << " is not in the valid range (ksPlus > 4.535) and now set to 0.0"<< std::endl; | ||
| 273 | ksPlus = 0.0; | ||
| 274 | printedRangeWarning = true; | ||
| 275 | } | ||
| 276 | ✗ | additionalRoughnessLength_[elementIdx] = 0.9 / (volVars.uStar() / volVars.kinematicViscosity()) | |
| 277 | ✗ | * (sqrt(ksPlus) - ksPlus * exp(-ksPlus / 6.0)); | |
| 278 | } | ||
| 279 | } | ||
| 280 | ✗ | } | |
| 281 | |||
| 282 | std::vector<Scalar> additionalRoughnessLength_; | ||
| 283 | std::vector<Scalar> kinematicEddyViscosity_; | ||
| 284 | |||
| 285 | //! Returns the implementation of the problem (i.e. static polymorphism) | ||
| 286 | Implementation &asImp_() | ||
| 287 | { return *static_cast<Implementation *>(this); } | ||
| 288 | |||
| 289 | //! \copydoc asImp_() | ||
| 290 | const Implementation &asImp_() const | ||
| 291 | { return *static_cast<const Implementation *>(this); } | ||
| 292 | }; | ||
| 293 | |||
| 294 | } // end namespace Dumux | ||
| 295 | |||
| 296 | #endif | ||
| 297 |