applications/utilities/preProcessing/wallFunctionTable/tabulatedWallFunction/SpaldingsLaw/SpaldingsLaw.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     namespace tabulatedWallFunctions
  34     {
  35         defineTypeNameAndDebug(SpaldingsLaw, 0);
  36         addToRunTimeSelectionTable
  37         (
  38             tabulatedWallFunction,
  39             SpaldingsLaw,
  40             dictionary
  41         );
  42     }
  43 }
  44
  45 const Foam::label Foam::tabulatedWallFunctions::SpaldingsLaw::maxIters_ = 1000;
  46
  47 const Foam::scalar
  48     Foam::tabulatedWallFunctions::SpaldingsLaw::tolerance_ = 1e-4;
  49
  50
  51 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  52
  53 void Foam::tabulatedWallFunctions::SpaldingsLaw::invertFunction()
  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         {
 107             WarningIn("SpaldingsLaw::invertFunction()")
 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
 120 Foam::tabulatedWallFunctions::SpaldingsLaw::SpaldingsLaw
 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
 141 Foam::tabulatedWallFunctions::SpaldingsLaw::~SpaldingsLaw()
 142 {}
 143
 144
 145 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 146
 147 Foam::scalar Foam::tabulatedWallFunctions::SpaldingsLaw::yPlus
 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
 160 Foam::scalar Foam::tabulatedWallFunctions::SpaldingsLaw::Re
 161 (
 162     const scalar uPlus
 163 ) const
 164 {
 165     return uPlus*yPlus(uPlus);
 166 }
 167
 168
 169 void Foam::tabulatedWallFunctions::SpaldingsLaw::writeData(Ostream& os) const
 170 {
 171     if (invertedTable_.log10())
 172     {
 173         os  << "log10(Re), y+, u+:" << endl;
 174         forAll(invertedTable_, i)
 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;
 185         forAll(invertedTable_, i)
 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 // ************************************************************************* //