src/turbulenceModels/RAS/compressible/wallFunctions/mutWallFunctions/mutStandardRoughWallFunction/mutStandardRoughWallFunctionFvPatchScalarField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "mutStandardRoughWallFunctionFvPatchScalarField.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 compressible
  38 {
  39 namespace RASModels
  40 {
  41
  42 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  43
  44 mutStandardRoughWallFunctionFvPatchScalarField::
  45 mutStandardRoughWallFunctionFvPatchScalarField
  46 (
  47     const fvPatch& p,
  48     const DimensionedField<scalar, volMesh>& iF
  49 )
  50 :
  51     fixedValueFvPatchScalarField(p, iF),
  52     roughnessHeight_(pTraits<scalar>::zero),
  53     roughnessConstant_(pTraits<scalar>::zero),
  54     roughnessFudgeFactor_(pTraits<scalar>::zero)
  55 {}
  56
  57
  58 mutStandardRoughWallFunctionFvPatchScalarField::
  59 mutStandardRoughWallFunctionFvPatchScalarField
  60 (
  61     const mutStandardRoughWallFunctionFvPatchScalarField& ptf,
  62     const fvPatch& p,
  63     const DimensionedField<scalar, volMesh>& iF,
  64     const fvPatchFieldMapper& mapper
  65 )
  66 :
  67     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
  68     roughnessHeight_(ptf.roughnessHeight_),
  69     roughnessConstant_(ptf.roughnessConstant_),
  70     roughnessFudgeFactor_(ptf.roughnessFudgeFactor_)
  71 {}
  72
  73
  74 mutStandardRoughWallFunctionFvPatchScalarField::
  75 mutStandardRoughWallFunctionFvPatchScalarField
  76 (
  77     const fvPatch& p,
  78     const DimensionedField<scalar, volMesh>& iF,
  79     const dictionary& dict
  80 )
  81 :
  82     fixedValueFvPatchScalarField(p, iF, dict),
  83     roughnessHeight_(readScalar(dict.lookup("roughnessHeight"))),
  84     roughnessConstant_(readScalar(dict.lookup("roughnessConstant"))),
  85     roughnessFudgeFactor_(readScalar(dict.lookup("roughnessFudgeFactor")))
  86 {}
  87
  88
  89 mutStandardRoughWallFunctionFvPatchScalarField::
  90 mutStandardRoughWallFunctionFvPatchScalarField
  91 (
  92     const mutStandardRoughWallFunctionFvPatchScalarField& tppsf
  93 )
  94 :
  95     fixedValueFvPatchScalarField(tppsf),
  96     roughnessHeight_(tppsf.roughnessHeight_),
  97     roughnessConstant_(tppsf.roughnessConstant_),
  98     roughnessFudgeFactor_(tppsf.roughnessFudgeFactor_)
  99 {}
 100
 101
 102 mutStandardRoughWallFunctionFvPatchScalarField::
 103 mutStandardRoughWallFunctionFvPatchScalarField
 104 (
 105     const mutStandardRoughWallFunctionFvPatchScalarField& tppsf,
 106     const DimensionedField<scalar, volMesh>& iF
 107 )
 108 :
 109     fixedValueFvPatchScalarField(tppsf, iF),
 110     roughnessHeight_(tppsf.roughnessHeight_),
 111     roughnessConstant_(tppsf.roughnessConstant_),
 112     roughnessFudgeFactor_(tppsf.roughnessFudgeFactor_)
 113 {}
 114
 115
 116 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 117
 118 void mutStandardRoughWallFunctionFvPatchScalarField::evaluate
 119 (
 120     const Pstream::commsTypes
 121 )
 122 {
 123     const RASModel& rasModel
 124         = db().lookupObject<RASModel>("RASProperties");
 125
 126     const scalar kappa = rasModel.kappa().value();
 127     const scalar E = rasModel.E().value();
 128     const scalar yPlusLam = 11.225;
 129
 130     // The reciprical of the distance to the adjacent cell centre.
 131     const scalarField& ry = patch().deltaCoeffs();
 132
 133     const fvPatchVectorField& U =
 134         patch().lookupPatchField<volVectorField, vector>("U");
 135
 136     const fvPatchScalarField& rho =
 137         patch().lookupPatchField<volScalarField, scalar>("rho");
 138
 139     // The flow velocity at the adjacent cell centre.
 140     scalarField magUp = mag(U.patchInternalField() - U);
 141
 142     const scalarField& muw =
 143         patch().lookupPatchField<volScalarField, scalar>("mu");
 144     scalarField& mutw = *this;
 145
 146     scalarField magFaceGradU = mag(U.snGrad());
 147
 148     if(roughnessHeight_ > 0.0)
 149     {
 150         // Rough Walls.
 151         const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_;
 152         static const scalar c_2 = 2.25/(90 - 2.25);
 153         static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25);
 154         static const scalar c_4 = c_3*log(2.25);
 155
 156         //if (KsPlusBasedOnYPlus_)
 157         {
 158             // If KsPlus is based on YPlus the extra term added to the law
 159             // of the wall will depend on yPlus.
 160             forAll(mutw, facei)
 161             {
 162                 const scalar magUpara = magUp[facei];
 163                 const scalar Re = rho[facei]*magUpara/(muw[facei]*ry[facei]);
 164                 const scalar kappaRe = kappa*Re;
 165
 166                 scalar yPlus = yPlusLam;
 167                 const scalar ryPlusLam = 1.0/yPlus;
 168
 169                 int iter = 0;
 170                 scalar yPlusLast = 0.0;
 171                 scalar dKsPlusdYPlus = roughnessHeight_*ry[facei];
 172
 173                 // Enforce the roughnessHeight to be less than the distance to
 174                 // the first cell centre.
 175                 if(dKsPlusdYPlus > 1)
 176                 {
 177                     dKsPlusdYPlus = 1;
 178                 }
 179
 180                 // Fudge factor to get results to be similar to fluent
 181                 // (at least difference between rough and smooth).
 182                 dKsPlusdYPlus *= roughnessFudgeFactor_;
 183
 184                 do
 185                 {
 186                     yPlusLast = yPlus;
 187
 188                     // The non-dimensional roughness height.
 189                     scalar KsPlus = yPlus*dKsPlusdYPlus;
 190
 191                     // The extra term in the law-of-the-wall.
 192                     scalar G = 0.0;
 193
 194                     scalar yPlusGPrime = 0.0;
 195
 196                     if (KsPlus >= 90)
 197                     {
 198                         const scalar t_1 = 1 + roughnessConstant_*KsPlus;
 199                         G = log(t_1);
 200                         yPlusGPrime = roughnessConstant_*KsPlus/t_1;
 201                     }
 202                     else if (KsPlus > 2.25)
 203                     {
 204                         const scalar t_1 = c_1*KsPlus - c_2;
 205                         const scalar t_2 = c_3*log(KsPlus) - c_4;
 206                         const scalar sint_2 = sin(t_2);
 207                         const scalar logt_1 = log(t_1);
 208                         G = logt_1*sint_2;
 209                         yPlusGPrime =
 210                             (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
 211                     }
 212
 213                     scalar denom = 1.0 + log(E* yPlus) - G - yPlusGPrime;
 214                     if(mag(denom) > VSMALL)
 215                     {
 216                         yPlus = (kappaRe + yPlus*(1 - yPlusGPrime))/denom;
 217                         if( yPlus < 0 )
 218                         {
 219                             yPlus = 0;
 220                         }
 221                     }
 222                     else
 223                     {
 224                         // Ensure immediate end and mutw = 0.
 225                         yPlus = 0;
 226                     }
 227
 228                 } while
 229                 (
 230                     mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001
 231                  && ++iter < 10
 232                  && yPlus > VSMALL
 233                 );
 234
 235                 if (yPlus > yPlusLam)
 236                 {
 237                     mutw[facei] = muw[facei]*(yPlus*yPlus/Re - 1);
 238                 }
 239                 else
 240                 {
 241                     mutw[facei] = 0.0;
 242                 }
 243             }
 244         }
 245     }
 246     else
 247     {
 248         // Smooth Walls.
 249         forAll(mutw, facei)
 250         {
 251             const scalar magUpara = magUp[facei];
 252             const scalar Re = rho[facei]*magUpara/(muw[facei]*ry[facei]);
 253             const scalar kappaRe = kappa*Re;
 254
 255             scalar yPlus = yPlusLam;
 256             const scalar ryPlusLam = 1.0/yPlus;
 257
 258             int iter = 0;
 259             scalar yPlusLast = 0.0;
 260
 261             do
 262             {
 263                 yPlusLast = yPlus;
 264                 yPlus = (kappaRe + yPlus)/(1.0 + log(E*yPlus));
 265
 266             } while
 267             (
 268                 mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001
 269              && ++iter < 10
 270             );
 271
 272             if (yPlus > yPlusLam)
 273             {
 274                 mutw[facei] = muw[facei]*(yPlus*yPlus/Re - 1);
 275             }
 276             else
 277             {
 278                 mutw[facei] = 0.0;
 279             }
 280         }
 281     }
 282 }
 283
 284
 285 void mutStandardRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
 286 {
 287     fixedValueFvPatchScalarField::write(os);
 288     os.writeKeyword("roughnessHeight")
 289         << roughnessHeight_ << token::END_STATEMENT << nl;
 290     os.writeKeyword("roughnessConstant")
 291         << roughnessConstant_ << token::END_STATEMENT << nl;
 292     os.writeKeyword("roughnessFudgeFactor")
 293         << roughnessFudgeFactor_ << token::END_STATEMENT << nl;
 294 }
 295
 296
 297 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 298
 299 makePatchTypeField
 300 (
 301     fvPatchScalarField,
 302     mutStandardRoughWallFunctionFvPatchScalarField
 303 );
 304
 305 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 306
 307 } // End namespace RASModels
 308 } // End namespace compressible
 309 } // End namespace Foam
 310
 311 // ************************************************************************* //