SpaldingsLaw.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2015 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "SpaldingsLaw.H"
30
31// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32
33namespace Foam
34{
35 namespace tabulatedWallFunctions
36 {
37 defineTypeNameAndDebug(SpaldingsLaw, 0);
39 (
40 tabulatedWallFunction,
41 SpaldingsLaw,
42 dictionary
43 );
44 }
45}
46
48
49const Foam::scalar
51
52
53// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
54
56{
57 // Initialise u+ and R
58 scalar Re = 0.0;
59 scalar uPlus = 1;
60
61 // Populate the table
63 {
65 {
66 Re = pow(10, (i*invertedTable_.dx() + invertedTable_.x0()));
67 }
68 else
69 {
71 }
72
73 // Use latest available u+ estimate
74 if (i > 0)
75 {
76 uPlus = invertedTable_[i-1];
77 }
78
79 // Newton iterations to determine u+
80 label iter = 0;
81 scalar error = GREAT;
82 do
83 {
84 scalar kUPlus = min(kappa_*uPlus, 50);
85 scalar A =
86 E_*sqr(uPlus)
87 + uPlus
88 *(exp(kUPlus) - pow3(kUPlus)/6 - 0.5*sqr(kUPlus) - kUPlus - 1);
89
90 scalar f = - Re + A/E_;
91
92 scalar df =
93 (
94 2*E_*uPlus
95 + exp(kUPlus)*(kUPlus + 1)
96 - 2.0/3.0*pow3(kUPlus)
97 - 1.5*sqr(kUPlus)
98 - 2*kUPlus
99 - 1
100 )/E_;
101
102 scalar uPlusNew = uPlus - f/(df + ROOTVSMALL);
103 error = mag((uPlus - uPlusNew)/uPlusNew);
104 uPlus = uPlusNew;
105 } while (error > tolerance_ && ++iter < maxIters_);
106
107 if (iter == maxIters_)
108 {
110 << "Newton iterations not converged:" << nl
111 << " iters = " << iter << ", error = " << error << endl;
112 }
113
114 // Set new values - constrain u+ >= 0
115 invertedTable_[i] = max(0, uPlus);
116 }
117}
118
119
120// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
121
123(
124 const dictionary& dict,
125 const polyMesh& mesh
126)
127:
128 tabulatedWallFunction(dict, mesh, typeName),
129 kappa_(coeffDict_.get<scalar>("kappa")),
130 E_(coeffDict_.get<scalar>("E"))
131{
132 invertFunction();
133
134 if (debug)
135 {
136 writeData(Info);
137 }
138}
139
140
141// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
142
144(
145 const scalar uPlus
146) const
147{
148 scalar kUPlus = min(kappa_*uPlus, 50);
149
150 return
151 uPlus
152 + 1/E_*(exp(kUPlus) - pow3(kUPlus)/6 - 0.5*sqr(kUPlus) - kUPlus - 1);
153}
154
155
157(
158 const scalar uPlus
159) const
160{
161 return uPlus*yPlus(uPlus);
162}
163
164
166{
167 if (invertedTable_.log10())
168 {
169 os << "log10(Re), y+, u+:" << endl;
170 forAll(invertedTable_, i)
171 {
172 scalar uPlus = invertedTable_[i];
173 scalar Re = ::log10(this->Re(uPlus));
174 scalar yPlus = this->yPlus(uPlus);
175 os << Re << ", " << yPlus << ", " << uPlus << endl;
176 }
177 }
178 else
179 {
180 os << "Re, y+, u+:" << endl;
181 forAll(invertedTable_, i)
182 {
183 scalar uPlus = invertedTable_[i];
184 scalar Re = this->Re(uPlus);
185 scalar yPlus = this->yPlus(uPlus);
186 os << Re << ", " << yPlus << ", " << uPlus << endl;
187 }
188 }
189}
190
191
192// ************************************************************************* //
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
scalar Re() const
Real part of complex number.
Definition: complexI.H:87
Computes the near wall for turbulence models.
Definition: yPlus.H:160
Computes U+ as a function of Reynolds number by inverting Spaldings law.
Definition: SpaldingsLaw.H:73
virtual void invertFunction()
Invert the function.
scalar kappa_
Von Karman constant.
Definition: SpaldingsLaw.H:79
static const scalar tolerance_
Tolerance.
Definition: SpaldingsLaw.H:91
static const label maxIters_
Maximum number of iterations.
Definition: SpaldingsLaw.H:88
scalar E_
Law-of-the-wall E coefficient.
Definition: SpaldingsLaw.H:82
virtual scalar Re(const scalar uPlus) const
Return Reynolds number as a function of u+.
virtual void writeData(Ostream &os) const
Write to Ostream.
uniformInterpolationTable< scalar > invertedTable_
Inverted table.
scalar dx() const
Return the fixed interval.
scalar x0() const
Return the lower limit.
const Switch & log10() const
Return the log10(x) flag.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
scalar uPlus
scalar yPlus
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar exp(const dimensionedScalar &ds)
scalarField Re(const UList< complex > &cf)
Extract real component.
Definition: complexField.C:159
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar log10(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
const bool writeData(pdfDictionary.get< bool >("writeData"))