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 | Copyright (C) 2011-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "SpaldingsLaw.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  namespace tabulatedWallFunctions
34  {
35  defineTypeNameAndDebug(SpaldingsLaw, 0);
37  (
38  tabulatedWallFunction,
39  SpaldingsLaw,
40  dictionary
41  );
42  }
43 }
44 
46 
47 const Foam::scalar
49 
50 
51 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
52 
54 {
55  // Initialise u+ and R
56  scalar Re = 0.0;
57  scalar uPlus = 1;
58 
59  // Populate the table
60  forAll(invertedTable_, i)
61  {
62  if (invertedTable_.log10())
63  {
64  Re = pow(10, (i*invertedTable_.dx() + invertedTable_.x0()));
65  }
66  else
67  {
68  Re = i*invertedTable_.dx() + invertedTable_.x0();
69  }
70 
71  // Use latest available u+ estimate
72  if (i > 0)
73  {
74  uPlus = invertedTable_[i-1];
75  }
76 
77  // Newton iterations to determine u+
78  label iter = 0;
79  scalar error = GREAT;
80  do
81  {
82  scalar kUPlus = min(kappa_*uPlus, 50);
83  scalar A =
84  E_*sqr(uPlus)
85  + uPlus
86  *(exp(kUPlus) - pow3(kUPlus)/6 - 0.5*sqr(kUPlus) - kUPlus - 1);
87 
88  scalar f = - Re + A/E_;
89 
90  scalar df =
91  (
92  2*E_*uPlus
93  + exp(kUPlus)*(kUPlus + 1)
94  - 2/3*pow3(kUPlus)
95  - 1.5*sqr(kUPlus)
96  - 2*kUPlus
97  - 1
98  )/E_;
99 
100  scalar uPlusNew = uPlus - f/(df + ROOTVSMALL);
101  error = mag((uPlus - uPlusNew)/uPlusNew);
102  uPlus = uPlusNew;
103  } while (error > tolerance_ && ++iter < maxIters_);
104 
105  if (iter == maxIters_)
106  {
108  << "Newton iterations not converged:" << nl
109  << " iters = " << iter << ", error = " << error << endl;
110  }
111 
112  // Set new values - constrain u+ >= 0
113  invertedTable_[i] = max(0, uPlus);
114  }
115 }
116 
117 
118 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
119 
121 (
122  const dictionary& dict,
123  const polyMesh& mesh
124 )
125 :
126  tabulatedWallFunction(dict, mesh, typeName),
127  kappa_(readScalar(coeffDict_.lookup("kappa"))),
128  E_(readScalar(coeffDict_.lookup("E")))
129 {
130  invertFunction();
131 
132  if (debug)
133  {
134  writeData(Info);
135  }
136 }
137 
138 
139 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
140 
142 {}
143 
144 
145 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
146 
148 (
149  const scalar uPlus
150 ) const
151 {
152  scalar kUPlus = min(kappa_*uPlus, 50);
153 
154  return
155  uPlus
156  + 1/E_*(exp(kUPlus) - pow3(kUPlus)/6 - 0.5*sqr(kUPlus) - kUPlus - 1);
157 }
158 
159 
161 (
162  const scalar uPlus
163 ) const
164 {
165  return uPlus*yPlus(uPlus);
166 }
167 
168 
170 {
171  if (invertedTable_.log10())
172  {
173  os << "log10(Re), y+, u+:" << endl;
175  {
176  scalar uPlus = invertedTable_[i];
177  scalar Re = ::log10(this->Re(uPlus));
178  scalar yPlus = this->yPlus(uPlus);
179  os << Re << ", " << yPlus << ", " << uPlus << endl;
180  }
181  }
182  else
183  {
184  os << "Re, y+, u+:" << endl;
186  {
187  scalar uPlus = invertedTable_[i];
188  scalar Re = this->Re(uPlus);
189  scalar yPlus = this->yPlus(uPlus);
190  os << Re << ", " << yPlus << ", " << uPlus << endl;
191  }
192  }
193 }
194 
195 
196 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const double e
Elementary charge.
Definition: doubleFloat.H:78
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
static const label maxIters_
Maximum number of iterations.
Definition: SpaldingsLaw.H:86
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual void invertFunction()
Invert the function.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Macros for easy insertion into run-time selection tables.
virtual void writeData(Ostream &os) const
Write to Ostream.
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedScalar exp(const dimensionedScalar &ds)
scalar E_
Law-of-the-wall E coefficient.
Definition: SpaldingsLaw.H:80
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
const bool writeData(readBool(pdfDictionary.lookup("writeData")))
const Switch & log10() const
Return the log10(x) flag.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual scalar Re(const scalar uPlus) const
Return Reynolds number as a function of u+.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
#define WarningInFunction
Report a warning using Foam::Warning.
virtual scalar yPlus(const scalar uPlus) const
Return y+ as a function of u+.
messageStream Info
static const scalar tolerance_
Tolerance.
Definition: SpaldingsLaw.H:89
dimensioned< scalar > mag(const dimensioned< Type > &)
uniformInterpolationTable< scalar > invertedTable_
Inverted table.
SpaldingsLaw(const dictionary &dict, const polyMesh &mesh)
dimensionedScalar log10(const dimensionedScalar &ds)
scalar kappa_
Von Karman constant.
Definition: SpaldingsLaw.H:77
Namespace for OpenFOAM.