src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H"
  28 #include "fvPatchFieldMapper.H"
  29 #include "volFields.H"
  30 #include "RASModel.H"
  31 #include "addToRunTimeSelectionTable.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37 namespace compressible
  38 {
  39 namespace RASModels
  40 {
  41
  42 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
  43
  44 tmp<scalarField>
  45 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
  46 (
  47     const scalarField& magUp
  48 ) const
  49 {
  50     const label patchI = patch().index();
  51
  52     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
  53     const scalarField& y = rasModel.y()[patchI];
  54     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
  55     const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
  56
  57     tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
  58     scalarField& yPlus = tyPlus();
  59
  60     if (roughnessHeight_ > 0.0)
  61     {
  62         // Rough Walls
  63         const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_;
  64         static const scalar c_2 = 2.25/(90 - 2.25);
  65         static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25);
  66         static const scalar c_4 = c_3*log(2.25);
  67
  68         //if (KsPlusBasedOnYPlus_)
  69         {
  70             // If KsPlus is based on YPlus the extra term added to the law
  71             // of the wall will depend on yPlus
  72             forAll(yPlus, facei)
  73             {
  74                 const scalar magUpara = magUp[facei];
  75                 const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
  76                 const scalar kappaRe = kappa_*Re;
  77
  78                 scalar yp = yPlusLam_;
  79                 const scalar ryPlusLam = 1.0/yp;
  80
  81                 int iter = 0;
  82                 scalar yPlusLast = 0.0;
  83                 scalar dKsPlusdYPlus = roughnessHeight_/y[facei];
  84
  85                 // Enforce the roughnessHeight to be less than the distance to
  86                 // the first cell centre.
  87                 if (dKsPlusdYPlus > 1)
  88                 {
  89                     dKsPlusdYPlus = 1;
  90                 }
  91
  92                 // Additional tuning parameter (fudge factor) - nominally = 1
  93                 dKsPlusdYPlus *= roughnessFudgeFactor_;
  94
  95                 do
  96                 {
  97                     yPlusLast = yp;
  98
  99                     // The non-dimensional roughness height
 100                     scalar KsPlus = yp*dKsPlusdYPlus;
 101
 102                     // The extra term in the law-of-the-wall
 103                     scalar G = 0.0;
 104
 105                     scalar yPlusGPrime = 0.0;
 106
 107                     if (KsPlus >= 90)
 108                     {
 109                         const scalar t_1 = 1 + roughnessConstant_*KsPlus;
 110                         G = log(t_1);
 111                         yPlusGPrime = roughnessConstant_*KsPlus/t_1;
 112                     }
 113                     else if (KsPlus > 2.25)
 114                     {
 115                         const scalar t_1 = c_1*KsPlus - c_2;
 116                         const scalar t_2 = c_3*log(KsPlus) - c_4;
 117                         const scalar sint_2 = sin(t_2);
 118                         const scalar logt_1 = log(t_1);
 119                         G = logt_1*sint_2;
 120                         yPlusGPrime =
 121                             (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
 122                     }
 123
 124                     scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime;
 125                     if (mag(denom) > VSMALL)
 126                     {
 127                         yp = (kappaRe + yp*(1 - yPlusGPrime))/denom;
 128                     }
 129
 130                 } while
 131                 (
 132                     mag(ryPlusLam*(yp - yPlusLast)) > 0.0001
 133                  && ++iter < 10
 134                  && yp > VSMALL
 135                 );
 136
 137                 yPlus[facei] = max(0.0, yp);
 138             }
 139         }
 140     }
 141     else
 142     {
 143         // Smooth Walls
 144         forAll(yPlus, facei)
 145         {
 146             const scalar magUpara = magUp[facei];
 147             const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
 148             const scalar kappaRe = kappa_*Re;
 149
 150             scalar yp = yPlusLam_;
 151             const scalar ryPlusLam = 1.0/yp;
 152
 153             int iter = 0;
 154             scalar yPlusLast = 0.0;
 155
 156             do
 157             {
 158                 yPlusLast = yp;
 159                 yp = (kappaRe + yp)/(1.0 + log(E_*yp));
 160
 161             } while
 162             (
 163                 mag(ryPlusLam*(yp - yPlusLast)) > 0.0001
 164              && ++iter < 10
 165             );
 166
 167             yPlus[facei] = max(0.0, yp);
 168         }
 169     }
 170
 171     return tyPlus;
 172 }
 173
 174
 175 tmp<scalarField>
 176 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
 177 {
 178     const label patchI = patch().index();
 179
 180     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
 181     const scalarField& y = rasModel.y()[patchI];
 182     const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
 183     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
 184     const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
 185
 186     scalarField magUp = mag(Uw.patchInternalField() - Uw);
 187
 188     tmp<scalarField> tyPlus = calcYPlus(magUp);
 189     scalarField& yPlus = tyPlus();
 190
 191     tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
 192     scalarField& mutw = tmutw();
 193
 194     forAll(yPlus, facei)
 195     {
 196         if (yPlus[facei] > yPlusLam_)
 197         {
 198             const scalar Re = rho[facei]*magUp[facei]*y[facei]/muw[facei];
 199             mutw[facei] = muw[facei]*(sqr(yPlus[facei])/Re - 1);
 200         }
 201     }
 202
 203     return tmutw;
 204 }
 205
 206
 207 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 208
 209 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
 210 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 211 (
 212     const fvPatch& p,
 213     const DimensionedField<scalar, volMesh>& iF
 214 )
 215 :
 216     mutWallFunctionFvPatchScalarField(p, iF),
 217     roughnessHeight_(pTraits<scalar>::zero),
 218     roughnessConstant_(pTraits<scalar>::zero),
 219     roughnessFudgeFactor_(pTraits<scalar>::zero)
 220 {}
 221
 222
 223 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
 224 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 225 (
 226     const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& ptf,
 227     const fvPatch& p,
 228     const DimensionedField<scalar, volMesh>& iF,
 229     const fvPatchFieldMapper& mapper
 230 )
 231 :
 232     mutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
 233     roughnessHeight_(ptf.roughnessHeight_),
 234     roughnessConstant_(ptf.roughnessConstant_),
 235     roughnessFudgeFactor_(ptf.roughnessFudgeFactor_)
 236 {}
 237
 238
 239 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
 240 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 241 (
 242     const fvPatch& p,
 243     const DimensionedField<scalar, volMesh>& iF,
 244     const dictionary& dict
 245 )
 246 :
 247     mutWallFunctionFvPatchScalarField(p, iF, dict),
 248     roughnessHeight_(readScalar(dict.lookup("roughnessHeight"))),
 249     roughnessConstant_(readScalar(dict.lookup("roughnessConstant"))),
 250     roughnessFudgeFactor_(readScalar(dict.lookup("roughnessFudgeFactor")))
 251 {}
 252
 253
 254 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
 255 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 256 (
 257     const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& rwfpsf
 258 )
 259 :
 260     mutWallFunctionFvPatchScalarField(rwfpsf),
 261     roughnessHeight_(rwfpsf.roughnessHeight_),
 262     roughnessConstant_(rwfpsf.roughnessConstant_),
 263     roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
 264 {}
 265
 266
 267 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
 268 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 269 (
 270     const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& rwfpsf,
 271     const DimensionedField<scalar, volMesh>& iF
 272 )
 273 :
 274     mutWallFunctionFvPatchScalarField(rwfpsf, iF),
 275     roughnessHeight_(rwfpsf.roughnessHeight_),
 276     roughnessConstant_(rwfpsf.roughnessConstant_),
 277     roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
 278 {}
 279
 280
 281 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 282
 283 tmp<scalarField>
 284 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
 285 {
 286     const label patchI = patch().index();
 287
 288     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
 289     const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
 290     const scalarField magUp = mag(Uw.patchInternalField() - Uw);
 291
 292     return calcYPlus(magUp);
 293 }
 294
 295
 296 void mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::write
 297 (
 298     Ostream& os
 299 ) const
 300 {
 301     fixedValueFvPatchScalarField::write(os);
 302     writeLocalEntries(os);
 303     os.writeKeyword("roughnessHeight")
 304         << roughnessHeight_ << token::END_STATEMENT << nl;
 305     os.writeKeyword("roughnessConstant")
 306         << roughnessConstant_ << token::END_STATEMENT << nl;
 307     os.writeKeyword("roughnessFudgeFactor")
 308         << roughnessFudgeFactor_ << token::END_STATEMENT << nl;
 309 }
 310
 311
 312 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 313
 314 makePatchTypeField
 315 (
 316     fvPatchScalarField,
 317     mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 318 );
 319
 320 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 321
 322 } // End namespace RASModels
 323 } // End namespace compressible
 324 } // End namespace Foam
 325
 326 // ************************************************************************* //