src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.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 "immersedBoundaryOmegaWallFunctionFvPatchScalarField.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 immersedBoundaryOmegaWallFunctionFvPatchScalarField::
  45 immersedBoundaryOmegaWallFunctionFvPatchScalarField
  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     beta1_(0.075)
  61 {}
  62
  63
  64 immersedBoundaryOmegaWallFunctionFvPatchScalarField::
  65 immersedBoundaryOmegaWallFunctionFvPatchScalarField
  66 (
  67     const fvPatch& p,
  68     const DimensionedField<scalar, volMesh>& iF,
  69     const dictionary& dict
  70 )
  71 :
  72     immersedBoundaryWallFunctionFvPatchScalarField(p, iF, dict),
  73     UName_(dict.lookupOrDefault<word>("U", "U")),
  74     kName_(dict.lookupOrDefault<word>("k", "k")),
  75     GName_(dict.lookupOrDefault<word>("G", "RASModel::G")),
  76     nuName_(dict.lookupOrDefault<word>("nu", "nu")),
  77     nutName_(dict.lookupOrDefault<word>("nut", "nut")),
  78     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
  79     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
  80     E_(dict.lookupOrDefault<scalar>("E", 9.8)),
  81     beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075))
  82 {}
  83
  84
  85 immersedBoundaryOmegaWallFunctionFvPatchScalarField::
  86 immersedBoundaryOmegaWallFunctionFvPatchScalarField
  87 (
  88     const immersedBoundaryOmegaWallFunctionFvPatchScalarField& ptf,
  89     const fvPatch& p,
  90     const DimensionedField<scalar, volMesh>& iF,
  91     const fvPatchFieldMapper& mapper
  92 )
  93 :
  94     immersedBoundaryWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
  95     UName_(ptf.UName_),
  96     kName_(ptf.kName_),
  97     GName_(ptf.GName_),
  98     nuName_(ptf.nuName_),
  99     nutName_(ptf.nutName_),
 100     Cmu_(ptf.Cmu_),
 101     kappa_(ptf.kappa_),
 102     E_(ptf.E_),
 103     beta1_(ptf.beta1_)
 104 {}
 105
 106
 107 immersedBoundaryOmegaWallFunctionFvPatchScalarField::
 108 immersedBoundaryOmegaWallFunctionFvPatchScalarField
 109 (
 110     const immersedBoundaryOmegaWallFunctionFvPatchScalarField& owfpsf
 111 )
 112 :
 113     immersedBoundaryWallFunctionFvPatchScalarField(owfpsf),
 114     UName_(owfpsf.UName_),
 115     kName_(owfpsf.kName_),
 116     GName_(owfpsf.GName_),
 117     nuName_(owfpsf.nuName_),
 118     nutName_(owfpsf.nutName_),
 119     Cmu_(owfpsf.Cmu_),
 120     kappa_(owfpsf.kappa_),
 121     E_(owfpsf.E_)
 122 {}
 123
 124
 125 immersedBoundaryOmegaWallFunctionFvPatchScalarField::
 126 immersedBoundaryOmegaWallFunctionFvPatchScalarField
 127 (
 128     const immersedBoundaryOmegaWallFunctionFvPatchScalarField& owfpsf,
 129     const DimensionedField<scalar, volMesh>& iF
 130 )
 131 :
 132     immersedBoundaryWallFunctionFvPatchScalarField(owfpsf, iF),
 133     UName_(owfpsf.UName_),
 134     kName_(owfpsf.kName_),
 135     GName_(owfpsf.GName_),
 136     nuName_(owfpsf.nuName_),
 137     nutName_(owfpsf.nutName_),
 138     Cmu_(owfpsf.Cmu_),
 139     kappa_(owfpsf.kappa_),
 140     E_(owfpsf.E_)
 141 {}
 142
 143
 144 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 145
 146 void immersedBoundaryOmegaWallFunctionFvPatchScalarField::updateCoeffs()
 147 {
 148     if (updated())
 149     {
 150         return;
 151     }
 152
 153     // If G field is not present, execute zero gradient evaluation
 154     // HJ, 20/Mar/2011
 155     if (!db().foundObject<volScalarField>(GName_))
 156     {
 157         InfoIn
 158         (
 159             "void immersedBoundaryOmegaWallFunctionFvPatchScalarField::"
 160             "updateCoeffs()"
 161         )   << "Cannot access " << GName_ << " field.  for patch "
 162             << patch().name() << ".  Evaluating as regular immersed boundary"
 163             << endl;
 164
 165         immersedBoundaryWallFunctionFvPatchScalarField::evaluate();
 166
 167         return;
 168     }
 169
 170     const vectorField& n = ibPatch().ibNormals();
 171
 172     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
 173     const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
 174
 175     const scalar Cmu25 = pow(Cmu_, 0.25);
 176     const scalar Cmu50 = sqrt(Cmu_);
 177
 178     volScalarField& G = const_cast<volScalarField&>
 179         (db().lookupObject<volScalarField>(GName_));
 180
 181     // Grab values of other fields required for wall functions
 182
 183     // Velocity
 184     const fvPatchVectorField& Uwg =
 185       patch().lookupPatchField<volVectorField, vector>(UName_);
 186     const immersedBoundaryVelocityWallFunctionFvPatchVectorField& Uw =
 187         refCast<const immersedBoundaryVelocityWallFunctionFvPatchVectorField>
 188         (
 189             Uwg
 190         );
 191
 192     // Calculate tangential component, taking into account wall velocity
 193     const scalarField UtanOld =
 194         mag((I - sqr(n)) & (Uw.ibSamplingPointValue() - Uw.ibValue()));
 195
 196     scalarField& UTangentialNew = Uw.wallTangentialValue();
 197
 198     // Turbulence kinetic energy
 199     const fvPatchScalarField& kg =
 200         patch().lookupPatchField<volScalarField, scalar>(kName_);
 201     const immersedBoundaryWallFunctionFvPatchScalarField& kw =
 202         refCast<const immersedBoundaryWallFunctionFvPatchScalarField>(kg);
 203
 204     // Current and new values of k at sampling point
 205     scalarField k = kw.ibSamplingPointValue();
 206     scalarField& kNew = kw.wallValue();
 207
 208     // Laminar viscosity
 209     const fvPatchScalarField& nuwg =
 210         patch().lookupPatchField<volScalarField, scalar>(nuName_);
 211     const immersedBoundaryFvPatchScalarField& nuw =
 212         refCast<const immersedBoundaryFvPatchScalarField>(nuwg);
 213     scalarField nu = nuw.ibCellValue();
 214
 215     // Turbulent viscosity
 216     const fvPatchScalarField& nutwg =
 217        patch().lookupPatchField<volScalarField, scalar>(nutName_);
 218     const immersedBoundaryWallFunctionFvPatchScalarField& nutw =
 219         refCast<const immersedBoundaryWallFunctionFvPatchScalarField>(nutwg);
 220
 221     // New values of nut
 222     scalarField nutOld = nutw.ibCellValue();
 223     scalarField& nutNew = nutw.wallValue();
 224
 225     const scalarField magGradUw = mag(Uw.ibGrad());
 226
 227     // Get the IB addressing and distance
 228     const labelList& ibc = ibPatch().ibCells();
 229
 230     // Distance to sampling point
 231     const scalarField& ySample = ibPatch().ibSamplingPointDelta();
 232
 233     // Distance from wall to IB point
 234     const scalarField& y = ibPatch().ibDelta();
 235
 236     // Omega: store IB cell values for direct insertion
 237     scalarField omegaSample = this->ibSamplingPointValue();
 238
 239     scalarField& omegaNew = this->wallValue();
 240
 241     // Mark values to be fixed
 242     boolList wf(ibc.size(), false);
 243
 244     // Calculate yPlus for sample points
 245     scalarField ypd = Cmu25*ySample*sqrt(k)/nu;
 246
 247     // Calculate wall function conditions
 248     forAll (ibc, ibCellI)
 249     {
 250         const scalar nuLam = nu[ibCellI];
 251
 252         // Calculate yPlus from k and laminar viscosity for the IB point
 253         const scalar yPlusSample = ypd[ibCellI];
 254
 255         scalar tauW, uTau;  // wall-shear and friction velocity from LOW
 256
 257         if (yPlusSample > yPlusLam)
 258         {
 259             // Calculate tauW from log-law using k and U at sampling point
 260
 261             tauW = UtanOld[ibCellI]*Cmu25*sqrt(k[ibCellI])*kappa_
 262                   /log(E_*yPlusSample);
 263         }
 264         else
 265         {
 266             // Sampling point is in laminar sublayer
 267             tauW = UtanOld[ibCellI]*Cmu25*sqrt(k[ibCellI])/yPlusSample;
 268         }
 269
 270         // friction velocity computed from k and U at sampling point
 271         uTau = sqrt(tauW);
 272
 273         // Calculate yPlus for IB point
 274
 275         scalar yPlusIB = yPlusSample*y[ibCellI]/ySample[ibCellI];
 276
 277         // Calculate wall function data in the immersed boundary point
 278         if (yPlusIB > yPlusLam)
 279         {
 280             // Logarithmic region
 281             wf[ibCellI] = true;
 282
 283             // turbulent viscosity at IB cell and at wall
 284             scalar nutw = nuLam*(yPlusIB*kappa_/log(E_*yPlusIB) - 1);
 285
 286             // Fix generation even though it if is not used
 287             G[ibc[ibCellI]] =
 288                 sqr((nutw + nuLam)*magGradUw[ibCellI])/
 289                 (Cmu25*sqrt(k[ibCellI])*kappa_*y[ibCellI]);
 290
 291             // Compute k at the IB cell
 292             kNew[ibCellI] = tauW/Cmu50;  // equilibrium boundary layer
 293             // kNew[ibCellI] = k[ibCellI];  // zero-Gradient (less stable)
 294
 295             // Compute omega at the IB cell
 296             omegaNew[ibCellI] = sqrt(kNew[ibCellI])/(Cmu25*kappa_*y[ibCellI]);
 297
 298             // Log-Law for tangential velocity  - uTau = Cmu25*sqrt(kNew)
 299             UTangentialNew[ibCellI] = uTau/kappa_*log(E_*yPlusIB);
 300
 301             // Calculate turbulent viscosity
 302             nutNew[ibCellI] = nutw;
 303         }
 304         else
 305         {
 306             // Laminar sub-layer
 307             wf[ibCellI] = false;
 308
 309             // G is zero  - immaterial!
 310             // G[ibc[ibCellI]] = 0;
 311
 312             // quadratic fit
 313             kNew[ibCellI] = k[ibCellI]*sqr(yPlusIB/yPlusLam);
 314
 315             // Compute omega at the IB cell
 316             omegaNew[ibCellI] = 6.0*nu[ibCellI]/(beta1_*sqr(y[ibCellI]));
 317
 318             // Laminar sub-layer for tangential velocity: uPlus = yPlus
 319             UTangentialNew[ibCellI] = uTau*yPlusIB;
 320
 321             // Turbulent viscosity is zero
 322             nutNew[ibCellI] = SMALL;
 323
 324         }
 325     }
 326
 327 //     Info<< "UTangentialNew " << min(UTangentialNew) << " " << max(UTangentialNew) << endl;
 328 //     Info<< "nutNew " << min(nutNew) << " " << max(nutNew) << endl;
 329 //     Info<< "kNew " << min(kNew) << " " << max(kNew) << endl;
 330 //     Info<< "epsilonNew " << min(epsilonNew) << " " << max(epsilonNew) << endl;
 331
 332     // Set the fields to calculated wall function values
 333     Uw.wallMask() = true;
 334     kw.wallMask() = wf;
 335     nutw.wallMask() = true;
 336     this->wallMask() = true;
 337
 338     // Insert epsilon values into the internal field
 339     immersedBoundaryWallFunctionFvPatchScalarField::updateCoeffs();
 340 }
 341
 342
 343 void immersedBoundaryOmegaWallFunctionFvPatchScalarField::evaluate
 344 (
 345     const Pstream::commsTypes commsType
 346 )
 347 {
 348     // Insert epsilon values into the internal field
 349     this->setIbCellValues(this->wallValue());
 350
 351     fvPatchScalarField::evaluate(commsType);
 352 }
 353
 354
 355 void immersedBoundaryOmegaWallFunctionFvPatchScalarField::
 356 write(Ostream& os) const
 357 {
 358     immersedBoundaryWallFunctionFvPatchScalarField::write(os);
 359     writeEntryIfDifferent<word>(os, "U", "U", UName_);
 360     writeEntryIfDifferent<word>(os, "k", "k", kName_);
 361     writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
 362     writeEntryIfDifferent<word>(os, "nu", "nu", nuName_);
 363     writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
 364     os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
 365     os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
 366     os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
 367     os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
 368 }
 369
 370
 371 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 372
 373 makePatchTypeField
 374 (
 375     fvPatchScalarField,
 376     immersedBoundaryOmegaWallFunctionFvPatchScalarField
 377 );
 378
 379 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 380
 381 } // End namespace RASModels
 382 } // End namespace incompressible
 383 } // End namespace Foam
 384
 385 // ************************************************************************* //