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