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