src/turbulenceModels/compressible/LES/derivedFvPatchFields/wallFunctions/alphaSgsWallFunctions/alphaSgsJayatillekeWallFunction/alphaSgsJayatillekeWallFunctionFvPatchScalarField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "alphaSgsJayatillekeWallFunctionFvPatchScalarField.H"
  28 #include "LESModel.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 LESModels
  40 {
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::maxExp_ = 50.0;
  45 scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::tolerance_ = 0.01;
  46 label alphaSgsJayatillekeWallFunctionFvPatchScalarField::maxIters_ = 10;
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 void alphaSgsJayatillekeWallFunctionFvPatchScalarField::checkType()
  51 {
  52     if (!patch().isWall())
  53     {
  54         FatalErrorIn
  55         (
  56             "alphaSgsJayatillekeWallFunctionFvPatchScalarField::checkType()"
  57         )
  58             << "Patch type for patch " << patch().name() << " must be wall\n"
  59             << "Current patch type is " << patch().type() << nl
  60             << exit(FatalError);
  61     }
  62 }
  63
  64
  65 scalar alphaSgsJayatillekeWallFunctionFvPatchScalarField::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 alphaSgsJayatillekeWallFunctionFvPatchScalarField::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 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
 109 alphaSgsJayatillekeWallFunctionFvPatchScalarField
 110 (
 111     const fvPatch& p,
 112     const DimensionedField<scalar, volMesh>& iF
 113 )
 114 :
 115     fixedValueFvPatchScalarField(p, iF),
 116     Prt_(0.85),
 117     kappa_(0.41),
 118     E_(9.8)
 119 {
 120     checkType();
 121 }
 122
 123
 124 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
 125 alphaSgsJayatillekeWallFunctionFvPatchScalarField
 126 (
 127     const alphaSgsJayatillekeWallFunctionFvPatchScalarField& ptf,
 128     const fvPatch& p,
 129     const DimensionedField<scalar, volMesh>& iF,
 130     const fvPatchFieldMapper& mapper
 131 )
 132 :
 133     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
 134     Prt_(ptf.Prt_),
 135     kappa_(ptf.kappa_),
 136     E_(ptf.E_)
 137 {}
 138
 139
 140 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
 141 alphaSgsJayatillekeWallFunctionFvPatchScalarField
 142 (
 143     const fvPatch& p,
 144     const DimensionedField<scalar, volMesh>& iF,
 145     const dictionary& dict
 146 )
 147 :
 148     fixedValueFvPatchScalarField(p, iF, dict),
 149     Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85)),
 150     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
 151     E_(dict.lookupOrDefault<scalar>("E", 9.8))
 152 {
 153     checkType();
 154 }
 155
 156
 157 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
 158 alphaSgsJayatillekeWallFunctionFvPatchScalarField
 159 (
 160     const alphaSgsJayatillekeWallFunctionFvPatchScalarField& awfpsf
 161 )
 162 :
 163     fixedValueFvPatchScalarField(awfpsf),
 164     Prt_(awfpsf.Prt_),
 165     kappa_(awfpsf.kappa_),
 166     E_(awfpsf.E_)
 167 {
 168     checkType();
 169 }
 170
 171
 172 alphaSgsJayatillekeWallFunctionFvPatchScalarField::
 173 alphaSgsJayatillekeWallFunctionFvPatchScalarField
 174 (
 175     const alphaSgsJayatillekeWallFunctionFvPatchScalarField& awfpsf,
 176     const DimensionedField<scalar, volMesh>& iF
 177 )
 178 :
 179     fixedValueFvPatchScalarField(awfpsf, iF),
 180     Prt_(awfpsf.Prt_),
 181     kappa_(awfpsf.kappa_),
 182     E_(awfpsf.E_)
 183 {
 184     checkType();
 185 }
 186
 187
 188 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 189
 190 void alphaSgsJayatillekeWallFunctionFvPatchScalarField::evaluate
 191 (
 192     const Pstream::commsTypes
 193 )
 194 {
 195     const LESModel& lesModel = db().lookupObject<LESModel>("LESProperties");
 196
 197     // Field data
 198     const label patchI = patch().index();
 199
 200     const scalarField& muw = lesModel.mu().boundaryField()[patchI];
 201     const scalarField muSgsw = lesModel.muSgs()().boundaryField()[patchI];
 202
 203     const scalarField& alphaw = lesModel.alpha().boundaryField()[patchI];
 204     scalarField& alphaSgsw = *this;
 205
 206     const fvPatchVectorField& Uw = lesModel.U().boundaryField()[patchI];
 207     const scalarField magUp = mag(Uw.patchInternalField() - Uw);
 208     const scalarField magGradUw = mag(Uw.snGrad());
 209
 210     const scalarField& rhow = lesModel.rho().boundaryField()[patchI];
 211     const fvPatchScalarField& hw =
 212         patch().lookupPatchField<volScalarField, scalar>("h");
 213
 214     const scalarField& ry = patch().deltaCoeffs();
 215
 216     // Heat flux [W/m2] - lagging alphaSgsw
 217     const scalarField qDot = (alphaw + alphaSgsw)*hw.snGrad();
 218
 219     // Populate boundary values
 220     forAll(alphaSgsw, faceI)
 221     {
 222         // Calculate uTau using Newton-Raphson iteration
 223         scalar uTau =
 224             sqrt((muSgsw[faceI] + muw[faceI])/rhow[faceI]*magGradUw[faceI]);
 225
 226         if (uTau > ROOTVSMALL)
 227         {
 228             label iter = 0;
 229             scalar err = GREAT;
 230
 231             do
 232             {
 233                 scalar kUu = min(kappa_*magUp[faceI]/uTau, maxExp_);
 234                 scalar fkUu = exp(kUu) - 1.0 - kUu*(1.0 + 0.5*kUu);
 235
 236                 scalar f =
 237                     - uTau/(ry[faceI]*muw[faceI]/rhow[faceI])
 238                     + magUp[faceI]/uTau
 239                     + 1.0/E_*(fkUu - 1.0/6.0*kUu*sqr(kUu));
 240
 241                 scalar df =
 242                     - 1.0/(ry[faceI]*muw[faceI]/rhow[faceI])
 243                     - magUp[faceI]/sqr(uTau)
 244                     - 1.0/E_*kUu*fkUu/uTau;
 245
 246                 scalar uTauNew = uTau - f/df;
 247                 err = mag((uTau - uTauNew)/uTau);
 248                 uTau = uTauNew;
 249
 250             } while (uTau>VSMALL && err>tolerance_ && ++iter<maxIters_);
 251
 252             scalar yPlus = uTau/ry[faceI]/(muw[faceI]/rhow[faceI]);
 253
 254             // Molecular Prandtl number
 255             scalar Pr = muw[faceI]/alphaw[faceI];
 256
 257             // Molecular-to-turbulenbt Prandtl number ratio
 258             scalar Prat = Pr/Prt_;
 259
 260             // Thermal sublayer thickness
 261             scalar P = Psmooth(Prat);
 262             scalar yPlusTherm = this->yPlusTherm(P, Prat);
 263
 264             // Evaluate new effective thermal diffusivity
 265             scalar alphaEff = 0.0;
 266             if (yPlus < yPlusTherm)
 267             {
 268                 scalar A = qDot[faceI]*rhow[faceI]*uTau/ry[faceI];
 269                 scalar B = qDot[faceI]*Pr*yPlus;
 270                 scalar C = Pr*0.5*rhow[faceI]*uTau*sqr(magUp[faceI]);
 271                 alphaEff = A/(B + C + VSMALL);
 272             }
 273             else
 274             {
 275                 scalar A = qDot[faceI]*rhow[faceI]*uTau/ry[faceI];
 276                 scalar B = qDot[faceI]*Prt_*(1.0/kappa_*log(E_*yPlus) + P);
 277                 scalar magUc = uTau/kappa_*log(E_*yPlusTherm) - mag(Uw[faceI]);
 278                 scalar C =
 279                     0.5*rhow[faceI]*uTau
 280                    *(Prt_*sqr(magUp[faceI]) + (Pr - Prt_)*sqr(magUc));
 281                 alphaEff = A/(B + C + VSMALL);
 282             }
 283
 284             // Update turbulent thermal diffusivity
 285             alphaSgsw[faceI] = max(0.0, alphaEff - alphaw[faceI]);
 286
 287             if (debug)
 288             {
 289                 Info<< "    uTau           = " << uTau << nl
 290                     << "    Pr             = " << Pr << nl
 291                     << "    Prt            = " << Prt_ << nl
 292                     << "    qDot           = " << qDot[faceI] << nl
 293                     << "    yPlus          = " << yPlus << nl
 294                     << "    yPlusTherm     = " << yPlusTherm << nl
 295                     << "    alphaEff       = " << alphaEff << nl
 296                     << "    alphaw         = " << alphaw[faceI] << nl
 297                     << "    alphaSgsw      = " << alphaSgsw[faceI] << nl
 298                     << endl;
 299             }
 300         }
 301         else
 302         {
 303             alphaSgsw[faceI] = 0.0;
 304         }
 305     }
 306 }
 307
 308
 309 void alphaSgsJayatillekeWallFunctionFvPatchScalarField::write(Ostream& os) const
 310 {
 311     fvPatchField<scalar>::write(os);
 312     os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
 313     os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
 314     os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
 315     writeEntry("value", os);
 316 }
 317
 318
 319 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 320
 321 makePatchTypeField
 322 (
 323     fvPatchScalarField,
 324     alphaSgsJayatillekeWallFunctionFvPatchScalarField
 325 );
 326
 327 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 328
 329 } // End namespace LESModels
 330 } // End namespace compressible
 331 } // End namespace Foam
 332
 333 // ************************************************************************* //