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