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