3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
freeflow/rans/zeroeq/problem.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_ZEROEQ_PROBLEM_HH
25#define DUMUX_ZEROEQ_PROBLEM_HH
26
27#include <string>
28
29#include <dune/common/math.hh>
37
38#include "model.hh"
39
40namespace Dumux {
41
49template<class TypeTag>
50class RANSProblemImpl<TypeTag, TurbulenceModel::zeroeq> : public RANSProblemBase<TypeTag>
51{
53 using Implementation = GetPropType<TypeTag, Properties::Problem>;
55 using Grid = typename GridView::Grid;
57
59 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
60 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
64 using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
66
67 enum {
68 dim = Grid::dimension,
69 };
70 using DimVector = Dune::FieldVector<Scalar, dim>;
71 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
72
73public:
79 RANSProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
80 : ParentType(gridGeometry, paramGroup)
81 {}
82
87 {
88 if (!ParentType::isFlatWallBounded())
89 {
90 DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, zero-eq models should only be used for flat channel geometries. "
91 << "\n If your geometry is a flat channel, please set the runtime parameter RANS.IsFlatWallBounded to true. \n");
92 }
93
94 ParentType::updateStaticWallProperties();
95
96 // update size and initial values of the global vectors
97 kinematicEddyViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
98 additionalRoughnessLength_.resize(this->gridGeometry().elementMapper().size(), 0.0);
99 }
100
108 void updateDynamicWallProperties(const SolutionVector& curSol)
109 {
110 ParentType::updateDynamicWallProperties(curSol);
111
112 // correct roughness lengths if a sand grain roughness is specified
113 if (hasParam("Problem.SandGrainRoughness"))
114 calculateRoughnessLength_(curSol);
115
116 // update routine for specfic models
117 if (eddyViscosityModel().compare("baldwinLomax") == 0)
118 updateBaldwinLomaxProperties();
119 }
120
125 {
126 std::vector<Scalar> kinematicEddyViscosityInner(this->gridGeometry().elementMapper().size(), 0.0);
127 std::vector<Scalar> kinematicEddyViscosityOuter(this->gridGeometry().elementMapper().size(), 0.0);
128 std::vector<Scalar> kinematicEddyViscosityDifference(this->gridGeometry().elementMapper().size(), 0.0);
129 std::vector<Scalar> switchingPosition(this->gridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max());
130
131 using std::abs;
132 using std::exp;
133 using std::min;
134 using std::sqrt;
135 using Dune::power;
136 const Scalar aPlus = 26.0;
137 const Scalar k = 0.0168;
138 const Scalar cCP = 1.6;
139 const Scalar cWake = 0.25;
140 const Scalar cKleb = 0.3;
141
142 std::vector<Scalar> storedFMax;
143 std::vector<Scalar> storedYFMax;
144 storedFMax.resize(this->gridGeometry().elementMapper().size(), 0.0);
145 storedYFMax.resize(this->gridGeometry().elementMapper().size(), 0.0);
146
147 // (1) calculate inner viscosity and Klebanoff function
148 for (const auto& element : elements(this->gridGeometry().gridView()))
149 {
150 unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
151 Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
152 unsigned int flowDirectionAxis = this->flowDirectionAxis(elementIdx);
153 unsigned int wallNormalAxis = this->wallNormalAxis(elementIdx);
154
155 Scalar omegaAbs = abs(this->velocityGradient(elementIdx, flowDirectionAxis, wallNormalAxis)
156 - this->velocityGradient(elementIdx, wallNormalAxis, flowDirectionAxis));
157 Scalar uStar = sqrt(this->kinematicViscosity(asImp_().wallElementIndex(elementIdx))
158 * abs(this->velocityGradient(asImp_().wallElementIndex(elementIdx), flowDirectionAxis, wallNormalAxis)));
159 Scalar yPlus = effectiveWallDistance * uStar / this->kinematicViscosity(elementIdx);
160 Scalar mixingLength = this->karmanConstant() * effectiveWallDistance * (1.0 - exp(-yPlus / aPlus));
161 kinematicEddyViscosityInner[elementIdx] = mixingLength * mixingLength * omegaAbs;
162
163 Scalar f = effectiveWallDistance * omegaAbs * (1.0 - exp(-yPlus / aPlus));
164 if (f > storedFMax[asImp_().wallElementIndex(elementIdx)])
165 {
166 storedFMax[asImp_().wallElementIndex(elementIdx)] = f;
167 storedYFMax[asImp_().wallElementIndex(elementIdx)] = effectiveWallDistance;
168 }
169 }
170
171 // (2) calculate outer viscosity
172 for (const auto& element : elements(this->gridGeometry().gridView()))
173 {
174 unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
175 Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
176
177 Scalar maxVelocityNorm = 0.0;
178 Scalar minVelocityNorm = 0.0;
179 for (unsigned axisIdx = 0; axisIdx < dim; ++axisIdx)
180 {
181 maxVelocityNorm += asImp_().velocityMaximum(asImp_().wallElementIndex(elementIdx))[axisIdx]
182 * asImp_().velocityMaximum(asImp_().wallElementIndex(elementIdx))[axisIdx];
183 minVelocityNorm += asImp_().velocityMinimum(asImp_().wallElementIndex(elementIdx))[axisIdx]
184 * asImp_().velocityMinimum(asImp_().wallElementIndex(elementIdx))[axisIdx];
185 }
186
187 Scalar deltaU = sqrt(maxVelocityNorm) - sqrt(minVelocityNorm);
188 Scalar yFMax = storedYFMax[asImp_().wallElementIndex(elementIdx)];
189 Scalar fMax = storedFMax[asImp_().wallElementIndex(elementIdx)];
190 Scalar fWake = min(yFMax * fMax, cWake * yFMax * deltaU * deltaU / fMax);
191 Scalar fKleb = 1.0 / (1.0 + 5.5 * power(cKleb * effectiveWallDistance / yFMax, 6));
192 kinematicEddyViscosityOuter[elementIdx] = k * cCP * fWake * fKleb;
193
194 kinematicEddyViscosityDifference[elementIdx]
195 = kinematicEddyViscosityInner[elementIdx] - kinematicEddyViscosityOuter[elementIdx];
196 }
197
198 // (3) switching point
199 for (const auto& element : elements(this->gridGeometry().gridView()))
200 {
201 unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
202 Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
203
204 // checks if sign switches, by multiplication
205 Scalar check = kinematicEddyViscosityDifference[asImp_().wallElementIndex(elementIdx)] * kinematicEddyViscosityDifference[elementIdx];
206 if (check < 0 // means sign has switched
207 && switchingPosition[asImp_().wallElementIndex(elementIdx)] > effectiveWallDistance)
208 {
209 switchingPosition[asImp_().wallElementIndex(elementIdx)] = effectiveWallDistance;
210 }
211 }
212
213 // (4) finally determine eddy viscosity
214 for (const auto& element : elements(this->gridGeometry().gridView()))
215 {
216 unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
217 Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength(elementIdx);
218
219 kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityInner[elementIdx];
220 if (effectiveWallDistance >= switchingPosition[asImp_().wallElementIndex(elementIdx)])
221 {
222 kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityOuter[elementIdx];
223 }
224 }
225 }
226
227 std::string eddyViscosityModel() const
228 {
229 static const std::string eddyViscosityModel = getParamFromGroup<std::string>(this->paramGroup(), "RANS.EddyViscosityModel", "vanDriest");
230 return eddyViscosityModel;
231 }
232
233 int additionalRoughnessLength(const int elementIdx) const
234 { return additionalRoughnessLength_[elementIdx]; }
235
236 Scalar kinematicEddyViscosity(const int elementIdx) const
237 { return kinematicEddyViscosity_[elementIdx]; }
238
239private:
240
241 void calculateRoughnessLength_(const SolutionVector& curSol)
242 {
243 bool printedRangeWarning = false;
244 auto fvGeometry = localView(this->gridGeometry());
245 for (const auto& element : elements(this->gridGeometry().gridView()))
246 {
247 static const Scalar sandGrainRoughness = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.SandGrainRoughness");
248 unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
249 fvGeometry.bindElement(element);
250
251 for (auto&& scv : scvs(fvGeometry))
252 {
253 using std::sqrt;
254 using std::exp;
255
256 const int dofIdx = scv.dofIndex();
257
258 // construct a privars object from the cell center solution vector
259 const auto& cellCenterPriVars = curSol[GridGeometry::cellCenterIdx()][dofIdx];
260 PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
261 auto elemSol = elementSolution<typename GridGeometry::LocalView>(std::move(priVars));
262
263 VolumeVariables volVars;
264 volVars.update(elemSol, asImp_(), element, scv);
265
266 Scalar ksPlus = sandGrainRoughness * volVars.uStar() / volVars.kinematicViscosity();
267 if (ksPlus > 0 && eddyViscosityModel().compare("baldwinLomax") == 0)
268 {
269 DUNE_THROW(Dune::NotImplemented, "Roughness is not implemented for the Baldwin-Lomax model.");
270 }
271 if (ksPlus > 2000.)
272 {
273 std::cout << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
274 << " is not in the valid range (ksPlus < 2000),"
275 << " for high ksPlus values the roughness function reaches a turning point."<< std::endl;
276 DUNE_THROW(Dune::InvalidStateException, "Unphysical roughness behavior.");
277 }
278 else if (ksPlus > 0.0 && ksPlus < 4.535 && !printedRangeWarning)
279 {
280 Dune::dinfo << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
281 << " is not in the valid range (ksPlus > 4.535) and now set to 0.0"<< std::endl;
282 ksPlus = 0.0;
283 printedRangeWarning = true;
284 }
285 additionalRoughnessLength_[elementIdx] = 0.9 / (volVars.uStar() / volVars.kinematicViscosity())
286 * (sqrt(ksPlus) - ksPlus * exp(-ksPlus / 6.0));
287 }
288 }
289 }
290
291 std::vector<Scalar> additionalRoughnessLength_;
292 std::vector<Scalar> kinematicEddyViscosity_;
293
295 Implementation &asImp_()
296 { return *static_cast<Implementation *>(this); }
297
299 const Implementation &asImp_() const
300 { return *static_cast<const Implementation *>(this); }
301};
302
303} // end namespace Dumux
304
305#endif
Base class for all staggered fv problems.
Free function to get the local view of a grid cache object.
The available discretization methods in Dumux.
The available free flow turbulence models in Dumux.
TurbulenceModel
The available free flow turbulence models in Dumux.
Definition: turbulencemodel.hh:38
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
bool hasParam(const std::string &param)
Check whether a key exists in the parameter tree.
Definition: parameters.hh:169
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
Definition: freeflow/navierstokes/problem.hh:55
forward declare
Definition: freeflow/rans/problem.hh:44
Reynolds-Averaged Navier-Stokes problem base class.
Definition: freeflow/rans/problem.hh:60
RANSProblemImpl(std::shared_ptr< const GridGeometry > gridGeometry, const std::string &paramGroup="")
The constructor.
Definition: freeflow/rans/zeroeq/problem.hh:79
void updateStaticWallProperties()
Correct size of the static (solution independent) wall variables.
Definition: freeflow/rans/zeroeq/problem.hh:86
int additionalRoughnessLength(const int elementIdx) const
Definition: freeflow/rans/zeroeq/problem.hh:233
void updateDynamicWallProperties(const SolutionVector &curSol)
Update the dynamic (solution dependent) relations to the walls.
Definition: freeflow/rans/zeroeq/problem.hh:108
Scalar kinematicEddyViscosity(const int elementIdx) const
Definition: freeflow/rans/zeroeq/problem.hh:236
std::string eddyViscosityModel() const
Definition: freeflow/rans/zeroeq/problem.hh:227
void updateBaldwinLomaxProperties()
Update the relations and coefficients for the Baldwin-Lomax turbulence model.
Definition: freeflow/rans/zeroeq/problem.hh:124
Declares all properties used in Dumux.
Adaption of the fully implicit scheme to the tracer transport model.
The local element solution class for staggered methods.