src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryEpsilonWallFunctions/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "immersedBoundaryEpsilonWallFunctionFvPatchScalarField.H"
  27 #include "immersedBoundaryVelocityWallFunctionFvPatchVectorField.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 incompressible
  38 {
  39 namespace RASModels
  40 {
  41
  42 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  43
  44 immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
  45 immersedBoundaryEpsilonWallFunctionFvPatchScalarField
  46 (
  47     const fvPatch& p,
  48     const DimensionedField<scalar, volMesh>& iF
  49 )
  50 :
  51     immersedBoundaryWallFunctionFvPatchScalarField(p, iF),
  52     UName_("U"),
  53     kName_("k"),
  54     GName_("RASModel::G"),
  55     nuName_("nu"),
  56     nutName_("nut"),
  57     Cmu_(0.09),
  58     kappa_(0.41),
  59     E_(9.8)
  60 {}
  61
  62
  63 immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
  64 immersedBoundaryEpsilonWallFunctionFvPatchScalarField
  65 (
  66     const fvPatch& p,
  67     const DimensionedField<scalar, volMesh>& iF,
  68     const dictionary& dict
  69 )
  70 :
  71     immersedBoundaryWallFunctionFvPatchScalarField(p, iF, dict),
  72     UName_(dict.lookupOrDefault<word>("U", "U")),
  73     kName_(dict.lookupOrDefault<word>("k", "k")),
  74     GName_(dict.lookupOrDefault<word>("G", "RASModel::G")),
  75     nuName_(dict.lookupOrDefault<word>("nu", "nu")),
  76     nutName_(dict.lookupOrDefault<word>("nut", "nut")),
  77     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
  78     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
  79     E_(dict.lookupOrDefault<scalar>("E", 9.8))
  80 {}
  81
  82
  83 immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
  84 immersedBoundaryEpsilonWallFunctionFvPatchScalarField
  85 (
  86     const immersedBoundaryEpsilonWallFunctionFvPatchScalarField& ptf,
  87     const fvPatch& p,
  88     const DimensionedField<scalar, volMesh>& iF,
  89     const fvPatchFieldMapper& mapper
  90 )
  91 :
  92     immersedBoundaryWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
  93     UName_(ptf.UName_),
  94     kName_(ptf.kName_),
  95     GName_(ptf.GName_),
  96     nuName_(ptf.nuName_),
  97     nutName_(ptf.nutName_),
  98     Cmu_(ptf.Cmu_),
  99     kappa_(ptf.kappa_),
 100     E_(ptf.E_)
 101 {}
 102
 103
 104 immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
 105 immersedBoundaryEpsilonWallFunctionFvPatchScalarField
 106 (
 107     const immersedBoundaryEpsilonWallFunctionFvPatchScalarField& ewfpsf
 108 )
 109 :
 110     immersedBoundaryWallFunctionFvPatchScalarField(ewfpsf),
 111     UName_(ewfpsf.UName_),
 112     kName_(ewfpsf.kName_),
 113     GName_(ewfpsf.GName_),
 114     nuName_(ewfpsf.nuName_),
 115     nutName_(ewfpsf.nutName_),
 116     Cmu_(ewfpsf.Cmu_),
 117     kappa_(ewfpsf.kappa_),
 118     E_(ewfpsf.E_)
 119 {}
 120
 121
 122 immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
 123 immersedBoundaryEpsilonWallFunctionFvPatchScalarField
 124 (
 125     const immersedBoundaryEpsilonWallFunctionFvPatchScalarField& ewfpsf,
 126     const DimensionedField<scalar, volMesh>& iF
 127 )
 128 :
 129     immersedBoundaryWallFunctionFvPatchScalarField(ewfpsf, iF),
 130     UName_(ewfpsf.UName_),
 131     kName_(ewfpsf.kName_),
 132     GName_(ewfpsf.GName_),
 133     nuName_(ewfpsf.nuName_),
 134     nutName_(ewfpsf.nutName_),
 135     Cmu_(ewfpsf.Cmu_),
 136     kappa_(ewfpsf.kappa_),
 137     E_(ewfpsf.E_)
 138 {}
 139
 140
 141 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 142
 143 void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::updateCoeffs()
 144 {
 145     if (updated())
 146     {
 147         return;
 148     }
 149
 150     // If G field is not present, execute zero gradient evaluation
 151     // HJ, 20/Mar/2011
 152     if (!db().foundObject<volScalarField>(GName_))
 153     {
 154         InfoIn
 155         (
 156             "void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::"
 157             "updateCoeffs()"
 158         )   << "Cannot access " << GName_ << " field.  for patch "
 159             << patch().name() << ".  Evaluating as regular immersed boundary"
 160             << endl;
 161
 162         immersedBoundaryWallFunctionFvPatchScalarField::evaluate();
 163
 164         return;
 165     }
 166
 167     const vectorField& n = ibPatch().ibNormals();
 168
 169     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
 170     const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
 171
 172     const scalar Cmu25 = pow(Cmu_, 0.25);
 173     const scalar Cmu50 = sqrt(Cmu_);
 174     const scalar Cmu75 = pow(Cmu_, 0.75);
 175
 176     volScalarField& G = const_cast<volScalarField&>
 177         (db().lookupObject<volScalarField>(GName_));
 178
 179     // Grab values of other fields required for wall functions
 180
 181     // Velocity
 182     const fvPatchVectorField& Uwg =
 183       patch().lookupPatchField<volVectorField, vector>(UName_);
 184     const immersedBoundaryVelocityWallFunctionFvPatchVectorField& Uw =
 185         refCast<const immersedBoundaryVelocityWallFunctionFvPatchVectorField>
 186         (
 187             Uwg
 188         );
 189
 190     // Calculate tangential component, taking into account wall velocity
 191     const vectorField UtanOld =
 192         (I - sqr(n)) & (Uw.ibSamplingPointValue() - Uw.ibValue());
 193     const scalarField magUtanOld = mag(UtanOld);
 194
 195     // Tangential velocity component
 196     scalarField& UTangentialNew = Uw.wallTangentialValue();
 197
 198     // Wall shear stress
 199     vectorField& tauWall = Uw.tauWall();
 200
 201     // Turbulence kinetic energy
 202     const fvPatchScalarField& kg =
 203         patch().lookupPatchField<volScalarField, scalar>(kName_);
 204     const immersedBoundaryWallFunctionFvPatchScalarField& kw =
 205         refCast<const immersedBoundaryWallFunctionFvPatchScalarField>(kg);
 206
 207     // Current and new values of k at sampling point
 208     scalarField k = kw.ibSamplingPointValue();
 209     scalarField& kNew = kw.wallValue();
 210
 211     // Laminar viscosity
 212     const fvPatchScalarField& nuwg =
 213         patch().lookupPatchField<volScalarField, scalar>(nuName_);
 214     const immersedBoundaryFvPatchScalarField& nuw =
 215         refCast<const immersedBoundaryFvPatchScalarField>(nuwg);
 216     scalarField nu = nuw.ibCellValue();
 217
 218     // Turbulent viscosity
 219     const fvPatchScalarField& nutwg =
 220        patch().lookupPatchField<volScalarField, scalar>(nutName_);
 221     const immersedBoundaryWallFunctionFvPatchScalarField& nutw =
 222         refCast<const immersedBoundaryWallFunctionFvPatchScalarField>(nutwg);
 223
 224     // New values of nut
 225     scalarField nutOld = nutw.ibCellValue();
 226     scalarField& nutNew = nutw.wallValue();
 227
 228     const scalarField magGradUw = mag(Uw.ibGrad());
 229
 230     // Get the IB addressing and distance
 231     const labelList& ibc = ibPatch().ibCells();
 232
 233     // Distance to sampling point
 234     const scalarField& ySample = ibPatch().ibSamplingPointDelta();
 235
 236     // Distance from wall to IB point
 237     const scalarField& y = ibPatch().ibDelta();
 238
 239     // Epsilon: store IB cell values for direct insertion
 240     scalarField epsilonSample = this->ibSamplingPointValue();
 241
 242     scalarField& epsilonNew = this->wallValue();
 243
 244     // Mark values to be fixed
 245     boolList wf(ibc.size(), false);
 246
 247     // Calculate yPlus for sample points
 248     scalarField ypd = Cmu25*ySample*sqrt(k)/nu;
 249
 250     // Calculate wall function conditions
 251     forAll (ibc, ibCellI)
 252     {
 253         const scalar nuLam = nu[ibCellI];
 254
 255         // Calculate yPlus from k and laminar viscosity for the IB point
 256         const scalar yPlusSample = ypd[ibCellI];
 257
 258         scalar uTau;
 259
 260         if (yPlusSample > yPlusLam)
 261         {
 262             // Calculate uTau from log-law, knowing sampled k and U
 263             uTau = magUtanOld[ibCellI]*kappa_/log(E_*yPlusSample);
 264         }
 265         else
 266         {
 267             // Sampling point is in laminar sublayer
 268             // Bug fix: HJ, 11/Aug/2014
 269             uTau = yPlusSample;
 270
 271         }
 272
 273         // Set wall shear stress
 274         tauWall[ibCellI] = sqr(uTau)*UtanOld[ibCellI]/(magUtanOld[ibCellI] + SMALL);
 275
 276         // Calculate yPlus for IB point
 277 //         scalar yPlusIB = uTau*y[ibCellI]/nuLam;
 278         scalar yPlusIB = yPlusSample*y[ibCellI]/ySample[ibCellI];
 279
 280         // Calculate wall function data in the immersed boundary point
 281         if (yPlusIB > yPlusLam)
 282         {
 283             // Logarithmic region
 284             wf[ibCellI] = true;
 285
 286             scalar nutw = nuLam*(yPlusIB*kappa_/log(E_*yPlusIB) - 1);
 287
 288             // Fix generation even though it if is not used
 289             G[ibc[ibCellI]] =
 290                 sqr((nutw + nuLam)*magGradUw[ibCellI])/
 291                 (Cmu25*sqrt(k[ibCellI])*kappa_*y[ibCellI]);
 292
 293             // Log-Law for tangential velocity
 294             UTangentialNew[ibCellI] =
 295                 min
 296                 (
 297                     magUtanOld[ibCellI],
 298                     uTau/kappa_*log(E_*yPlusIB)
 299                 );
 300
 301             // Calculate turbulent viscosity
 302             nutNew[ibCellI] = nutw;
 303
 304             // Calculate k in the IB cell from G = epsilon
 305             kNew[ibCellI] = (nutw + nuLam)*magGradUw[ibCellI]/Cmu50;
 306
 307             // Calculate epsilon from yPlus and set it
 308             epsilonNew[ibCellI] =
 309                 Cmu75*pow(kNew[ibCellI], 1.5)/(kappa_*y[ibCellI]);
 310         }
 311         else
 312         {
 313             // Laminar sub-layer
 314             wf[ibCellI] = false;
 315
 316             // G is zero
 317             G[ibc[ibCellI]] = 0;
 318
 319             // Laminar sub-layer for tangential velocity: uPlus = yPlus
 320             UTangentialNew[ibCellI] = min(magUtanOld[ibCellI], uTau*yPlusIB);
 321
 322             // Turbulent viscosity is zero
 323             nutNew[ibCellI] = SMALL;
 324
 325             // k is zero gradient: use the sampled value
 326             kNew[ibCellI] = k[ibCellI];
 327
 328             // Calculate epsilon from yPlus and set it.
 329             // Note:  calculating equilibrium epsilon in the sub-layer creates
 330             //        an unrealistic oscillation: use sampled value
 331             // HJ, 27/Jul/2012
 332             epsilonNew[ibCellI] = epsilonSample[ibCellI];
 333         }
 334     }
 335
 336 //     Info<< "UTangentialNew " << min(UTangentialNew) << " " << max(UTangentialNew) << endl;
 337 //     Info<< "nutNew " << min(nutNew) << " " << max(nutNew) << endl;
 338 //     Info<< "kNew " << min(kNew) << " " << max(kNew) << endl;
 339 //     Info<< "epsilonNew " << min(epsilonNew) << " " << max(epsilonNew) << endl;
 340
 341     // Set the fields to calculated wall function values
 342     Uw.wallMask() = true;
 343     kw.wallMask() = wf;
 344     nutw.wallMask() = true;
 345     this->wallMask() = true;
 346
 347     // Insert epsilon values into the internal field
 348     immersedBoundaryWallFunctionFvPatchScalarField::updateCoeffs();
 349 }
 350
 351
 352 void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::evaluate
 353 (
 354     const Pstream::commsTypes commsType
 355 )
 356 {
 357     // Insert epsilon values into the internal field
 358     this->setIbCellValues(this->wallValue());
 359
 360     fvPatchScalarField::evaluate(commsType);
 361 }
 362
 363
 364 void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
 365 write(Ostream& os) const
 366 {
 367     immersedBoundaryWallFunctionFvPatchScalarField::write(os);
 368     writeEntryIfDifferent<word>(os, "U", "U", UName_);
 369     writeEntryIfDifferent<word>(os, "k", "k", kName_);
 370     writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
 371     writeEntryIfDifferent<word>(os, "nu", "nu", nuName_);
 372     writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
 373     os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
 374     os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
 375     os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
 376 }
 377
 378
 379 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 380
 381 makePatchTypeField
 382 (
 383     fvPatchScalarField,
 384     immersedBoundaryEpsilonWallFunctionFvPatchScalarField
 385 );
 386
 387 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 388
 389 } // End namespace RASModels
 390 } // End namespace incompressible
 391 } // End namespace Foam
 392
 393 // ************************************************************************* //