src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
   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 "kappatJayatillekeWallFunctionFvPatchScalarField.H"
  27 #include "RASModel.H"
  28 #include "fvPatchFieldMapper.H"
  29 #include "volFields.H"
  30 #include "wallFvPatch.H"
  31 #include "addToRunTimeSelectionTable.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37 namespace incompressible
  38 {
  39 namespace RASModels
  40 {
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 scalar kappatJayatillekeWallFunctionFvPatchScalarField::tolerance_ = 0.01;
  45 label kappatJayatillekeWallFunctionFvPatchScalarField::maxIters_ = 10;
  46
  47 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
  48
  49 void kappatJayatillekeWallFunctionFvPatchScalarField::checkType()
  50 {
  51     if (!isA<wallFvPatch>(patch()))
  52     {
  53         FatalErrorIn
  54         (
  55             "kappatJayatillekeWallFunctionFvPatchScalarField::checkType()"
  56         )   << "Invalid wall function specification" << nl
  57             << "    Patch type for patch " << patch().name()
  58             << " must be wall" << nl
  59             << "    Current patch type is " << patch().type() << nl << endl
  60             << abort(FatalError);
  61     }
  62 }
  63
  64
  65 scalar kappatJayatillekeWallFunctionFvPatchScalarField::Psmooth
  66 (
  67     const scalar Prat
  68 ) const
  69 {
  70     return 9.24*(pow(Prat, 0.75) - 1.0)*(1.0 + 0.28*exp(-0.007*Prat));
  71 }
  72
  73
  74 scalar kappatJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
  75 (
  76     const scalar P,
  77     const scalar Prat
  78 ) const
  79 {
  80     scalar ypt = 11.0;
  81
  82     for (int i=0; i<maxIters_; i++)
  83     {
  84         scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat;
  85         scalar df = 1.0 - 1.0/(ypt*kappa_*Prat);
  86         scalar yptNew = ypt - f/df;
  87
  88         if (yptNew < VSMALL)
  89         {
  90             return 0;
  91         }
  92         else if (mag(yptNew - ypt) < tolerance_)
  93         {
  94             return yptNew;
  95         }
  96         else
  97         {
  98             ypt = yptNew;
  99         }
 100      }
 101
 102     return ypt;
 103 }
 104
 105
 106 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 107
 108 kappatJayatillekeWallFunctionFvPatchScalarField::
 109 kappatJayatillekeWallFunctionFvPatchScalarField
 110 (
 111     const fvPatch& p,
 112     const DimensionedField<scalar, volMesh>& iF
 113 )
 114 :
 115     fixedValueFvPatchScalarField(p, iF),
 116     Prt_(0.85),
 117     Cmu_(0.09),
 118     kappa_(0.41),
 119     E_(9.8)
 120 {
 121     checkType();
 122 }
 123
 124
 125 kappatJayatillekeWallFunctionFvPatchScalarField::
 126 kappatJayatillekeWallFunctionFvPatchScalarField
 127 (
 128     const kappatJayatillekeWallFunctionFvPatchScalarField& ptf,
 129     const fvPatch& p,
 130     const DimensionedField<scalar, volMesh>& iF,
 131     const fvPatchFieldMapper& mapper
 132 )
 133 :
 134     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
 135     Prt_(ptf.Prt_),
 136     Cmu_(ptf.Cmu_),
 137     kappa_(ptf.kappa_),
 138     E_(ptf.E_)
 139 {
 140     checkType();
 141 }
 142
 143
 144 kappatJayatillekeWallFunctionFvPatchScalarField::
 145 kappatJayatillekeWallFunctionFvPatchScalarField
 146 (
 147     const fvPatch& p,
 148     const DimensionedField<scalar, volMesh>& iF,
 149     const dictionary& dict
 150 )
 151 :
 152     fixedValueFvPatchScalarField(p, iF, dict),
 153     Prt_(readScalar(dict.lookup("Prt"))), // force read to avoid ambiguity
 154     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
 155     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
 156     E_(dict.lookupOrDefault<scalar>("E", 9.8))
 157 {
 158     checkType();
 159 }
 160
 161
 162 kappatJayatillekeWallFunctionFvPatchScalarField::
 163 kappatJayatillekeWallFunctionFvPatchScalarField
 164 (
 165     const kappatJayatillekeWallFunctionFvPatchScalarField& wfpsf
 166 )
 167 :
 168     fixedValueFvPatchScalarField(wfpsf),
 169     Prt_(wfpsf.Prt_),
 170     Cmu_(wfpsf.Cmu_),
 171     kappa_(wfpsf.kappa_),
 172     E_(wfpsf.E_)
 173 {
 174     checkType();
 175 }
 176
 177
 178 kappatJayatillekeWallFunctionFvPatchScalarField::
 179 kappatJayatillekeWallFunctionFvPatchScalarField
 180 (
 181     const kappatJayatillekeWallFunctionFvPatchScalarField& wfpsf,
 182     const DimensionedField<scalar, volMesh>& iF
 183 )
 184 :
 185     fixedValueFvPatchScalarField(wfpsf, iF),
 186     Prt_(wfpsf.Prt_),
 187     Cmu_(wfpsf.Cmu_),
 188     kappa_(wfpsf.kappa_),
 189     E_(wfpsf.E_)
 190 {
 191     checkType();
 192 }
 193
 194
 195 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 196
 197 void kappatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
 198 {
 199     if (updated())
 200     {
 201         return;
 202     }
 203
 204     const label patchI = patch().index();
 205
 206     // Retrieve turbulence properties from model
 207     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
 208     const scalar Cmu25 = pow(Cmu_, 0.25);
 209     const scalarField& y = rasModel.y()[patchI];
 210     const scalarField& nuw = rasModel.nu()().boundaryField()[patchI];
 211     const tmp<volScalarField> tk = rasModel.k();
 212     const volScalarField& k = tk();
 213
 214     // Molecular Prandtl number
 215     const scalar
 216         Pr(dimensionedScalar(rasModel.transport().lookup("Pr")).value());
 217
 218     // Populate boundary values
 219     scalarField& kappatw = *this;
 220     forAll(kappatw, faceI)
 221     {
 222         label faceCellI = patch().faceCells()[faceI];
 223
 224         // y+
 225         scalar yPlus = Cmu25*sqrt(k[faceCellI])*y[faceI]/nuw[faceI];
 226
 227         // Molecular-to-turbulent Prandtl number ratio
 228         scalar Prat = Pr/Prt_;
 229
 230         // Thermal sublayer thickness
 231         scalar P = Psmooth(Prat);
 232         scalar yPlusTherm = this->yPlusTherm(P, Prat);
 233
 234         // Update turbulent thermal conductivity
 235         if (yPlus > yPlusTherm)
 236         {
 237             scalar nu = nuw[faceI];
 238             scalar kt = nu*(yPlus/(Prt_*(log(E_*yPlus)/kappa_ + P)) - 1/Pr);
 239             kappatw[faceI] = max(0.0, kt);
 240         }
 241         else
 242         {
 243             kappatw[faceI] = 0.0;
 244         }
 245     }
 246
 247     fixedValueFvPatchField<scalar>::updateCoeffs();
 248 }
 249
 250
 251 void kappatJayatillekeWallFunctionFvPatchScalarField::write(Ostream& os) const
 252 {
 253     fvPatchField<scalar>::write(os);
 254     os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
 255     os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
 256     os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
 257     os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
 258     writeEntry("value", os);
 259 }
 260
 261
 262 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 263
 264 makePatchTypeField
 265 (
 266     fvPatchScalarField,
 267     kappatJayatillekeWallFunctionFvPatchScalarField
 268 );
 269
 270 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 271
 272 } // End namespace RASModels
 273 } // End namespace incompressible
 274 } // End namespace Foam
 275
 276 // ************************************************************************* //