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