src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.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 "kOmegaSST.H"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 #include "backwardsCompatibilityWallFunctions.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35 namespace incompressible
  36 {
  37 namespace RASModels
  38 {
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(kOmegaSST, 0);
  43 addToRunTimeSelectionTable(RASModel, kOmegaSST, dictionary);
  44
  45 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  46
  47 tmp<volScalarField> kOmegaSST::F1(const volScalarField& CDkOmega) const
  48 {
  49     volScalarField CDkOmegaPlus = max
  50     (
  51         CDkOmega,
  52         dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
  53     );
  54
  55     volScalarField arg1 = min
  56     (
  57         min
  58         (
  59             max
  60             (
  61                 (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
  62                 scalar(500)*nu()/(sqr(y_)*omega_)
  63             ),
  64             (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
  65         ),
  66         scalar(10)
  67     );
  68
  69     return tanh(pow4(arg1));
  70 }
  71
  72
  73 tmp<volScalarField> kOmegaSST::F2() const
  74 {
  75     volScalarField arg2 = min
  76     (
  77         max
  78         (
  79             (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
  80             scalar(500)*nu()/(sqr(y_)*omega_)
  81         ),
  82         scalar(100)
  83     );
  84
  85     return tanh(sqr(arg2));
  86 }
  87
  88
  89 tmp<volScalarField> kOmegaSST::F3() const
  90 {
  91     tmp<volScalarField> arg3 = min
  92     (
  93         150*nu()/(omega_*sqr(y_)),
  94         scalar(10)
  95     );
  96
  97     return 1 - tanh(pow4(arg3));
  98 }
  99
 100
 101 tmp<volScalarField> kOmegaSST::F23() const
 102 {
 103     tmp<volScalarField> f23(F2());
 104
 105     if (F3_)
 106     {
 107         f23() *= F3();
 108     }
 109
 110     return f23;
 111 }
 112
 113
 114 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 115
 116 kOmegaSST::kOmegaSST
 117 (
 118     const volVectorField& U,
 119     const surfaceScalarField& phi,
 120     transportModel& lamTransportModel
 121 )
 122 :
 123     RASModel(typeName, U, phi, lamTransportModel),
 124
 125     alphaK1_
 126     (
 127         dimensionedScalar::lookupOrAddToDict
 128         (
 129             "alphaK1",
 130             coeffDict_,
 131             0.85
 132         )
 133     ),
 134     alphaK2_
 135     (
 136         dimensionedScalar::lookupOrAddToDict
 137         (
 138             "alphaK2",
 139             coeffDict_,
 140             1.0
 141         )
 142     ),
 143     alphaOmega1_
 144     (
 145         dimensionedScalar::lookupOrAddToDict
 146         (
 147             "alphaOmega1",
 148             coeffDict_,
 149             0.5
 150         )
 151     ),
 152     alphaOmega2_
 153     (
 154         dimensionedScalar::lookupOrAddToDict
 155         (
 156             "alphaOmega2",
 157             coeffDict_,
 158             0.856
 159         )
 160     ),
 161     gamma1_
 162     (
 163         dimensionedScalar::lookupOrAddToDict
 164         (
 165             "gamma1",
 166             coeffDict_,
 167             5.0/9.0
 168         )
 169     ),
 170     gamma2_
 171     (
 172         dimensionedScalar::lookupOrAddToDict
 173         (
 174             "gamma2",
 175             coeffDict_,
 176             0.44
 177         )
 178     ),
 179     beta1_
 180     (
 181         dimensionedScalar::lookupOrAddToDict
 182         (
 183             "beta1",
 184             coeffDict_,
 185             0.075
 186         )
 187     ),
 188     beta2_
 189     (
 190         dimensionedScalar::lookupOrAddToDict
 191         (
 192             "beta2",
 193             coeffDict_,
 194             0.0828
 195         )
 196     ),
 197     betaStar_
 198     (
 199         dimensionedScalar::lookupOrAddToDict
 200         (
 201             "betaStar",
 202             coeffDict_,
 203             0.09
 204         )
 205     ),
 206     a1_
 207     (
 208         dimensionedScalar::lookupOrAddToDict
 209         (
 210             "a1",
 211             coeffDict_,
 212             0.31
 213         )
 214     ),
 215     b1_
 216     (
 217         dimensionedScalar::lookupOrAddToDict
 218         (
 219             "b1",
 220             coeffDict_,
 221             1.0
 222         )
 223     ),
 224     c1_
 225     (
 226         dimensionedScalar::lookupOrAddToDict
 227         (
 228             "c1",
 229             coeffDict_,
 230             10.0
 231         )
 232     ),
 233     F3_
 234     (
 235         Switch::lookupOrAddToDict
 236         (
 237             "F3",
 238             coeffDict_,
 239             false
 240         )
 241     ),
 242
 243     y_(mesh_),
 244
 245     k_
 246     (
 247         IOobject
 248         (
 249             "k",
 250             runTime_.timeName(),
 251             U_.db(),
 252             IOobject::NO_READ,
 253             IOobject::AUTO_WRITE
 254         ),
 255         autoCreateK("k", mesh_, U_.db())
 256     ),
 257     omega_
 258     (
 259         IOobject
 260         (
 261             "omega",
 262             runTime_.timeName(),
 263             U_.db(),
 264             IOobject::NO_READ,
 265             IOobject::AUTO_WRITE
 266         ),
 267         autoCreateOmega("omega", mesh_, U_.db())
 268     ),
 269     nut_
 270     (
 271         IOobject
 272         (
 273             "nut",
 274             runTime_.timeName(),
 275             U_.db(),
 276             IOobject::NO_READ,
 277             IOobject::AUTO_WRITE
 278         ),
 279         autoCreateNut("nut", mesh_, U_.db())
 280     )
 281 {
 282     bound(k_, k0_);
 283     bound(omega_, omega0_);
 284
 285     nut_ =
 286     (
 287         a1_*k_/
 288         max
 289         (
 290             a1_*omega_,
 291             b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
 292         )
 293     );
 294     nut_.correctBoundaryConditions();
 295
 296     printCoeffs();
 297 }
 298
 299
 300 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 301
 302 tmp<volSymmTensorField> kOmegaSST::R() const
 303 {
 304     return tmp<volSymmTensorField>
 305     (
 306         new volSymmTensorField
 307         (
 308             IOobject
 309             (
 310                 "R",
 311                 runTime_.timeName(),
 312                 U_.db(),
 313                 IOobject::NO_READ,
 314                 IOobject::NO_WRITE
 315             ),
 316             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 317             k_.boundaryField().types()
 318         )
 319     );
 320 }
 321
 322
 323 tmp<volSymmTensorField> kOmegaSST::devReff() const
 324 {
 325     return tmp<volSymmTensorField>
 326     (
 327         new volSymmTensorField
 328         (
 329             IOobject
 330             (
 331                 "devRhoReff",
 332                 runTime_.timeName(),
 333                 U_.db(),
 334                 IOobject::NO_READ,
 335                 IOobject::NO_WRITE
 336             ),
 337            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 338         )
 339     );
 340 }
 341
 342
 343 tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
 344 {
 345     return
 346     (
 347       - fvm::laplacian(nuEff(), U)
 348       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 349     );
 350 }
 351
 352
 353 bool kOmegaSST::read()
 354 {
 355     if (RASModel::read())
 356     {
 357         alphaK1_.readIfPresent(coeffDict());
 358         alphaK2_.readIfPresent(coeffDict());
 359         alphaOmega1_.readIfPresent(coeffDict());
 360         alphaOmega2_.readIfPresent(coeffDict());
 361         gamma1_.readIfPresent(coeffDict());
 362         gamma2_.readIfPresent(coeffDict());
 363         beta1_.readIfPresent(coeffDict());
 364         beta2_.readIfPresent(coeffDict());
 365         betaStar_.readIfPresent(coeffDict());
 366         a1_.readIfPresent(coeffDict());
 367         b1_.readIfPresent(coeffDict());
 368         c1_.readIfPresent(coeffDict());
 369         F3_.readIfPresent("F3", coeffDict());
 370
 371         return true;
 372     }
 373     else
 374     {
 375         return false;
 376     }
 377 }
 378
 379
 380 void kOmegaSST::correct()
 381 {
 382     // Bound in case of topological change
 383     // HJ, 22/Aug/2007
 384     if (mesh_.changing())
 385     {
 386         bound(k_, k0_);
 387         bound(omega_, omega0_);
 388     }
 389
 390     RASModel::correct();
 391
 392     if (!turbulence_)
 393     {
 394         return;
 395     }
 396
 397     if (mesh_.changing())
 398     {
 399         y_.correct();
 400     }
 401
 402     const volScalarField S2(2*magSqr(symm(fvc::grad(U_))));
 403     volScalarField G("RASModel::G", nut_*S2);
 404
 405     // Update omega and G at the wall
 406     omega_.boundaryField().updateCoeffs();
 407
 408     const volScalarField CDkOmega
 409     (
 410         (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
 411     );
 412
 413     const volScalarField F1(this->F1(CDkOmega));
 414
 415     // Turbulent frequency equation
 416     fvScalarMatrix omegaEqn
 417     (
 418         fvm::ddt(omega_)
 419       + fvm::div(phi_, omega_)
 420       + fvm::SuSp(-fvc::div(phi_), omega_)
 421       - fvm::laplacian(DomegaEff(F1), omega_)
 422      ==
 423         gamma(F1)
 424        *min
 425         (
 426             S2,
 427             (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
 428         )
 429       - fvm::Sp(beta(F1)*omega_, omega_)
 430       - fvm::SuSp
 431         (
 432             (F1 - scalar(1))*CDkOmega/omega_,
 433             omega_
 434         )
 435     );
 436
 437     omegaEqn.relax();
 438
 439     // No longer needed: matrix completes at the point of solution
 440     // HJ, 17/Apr/2012
 441 //     omegaEqn.completeAssembly();
 442
 443     solve(omegaEqn);
 444     bound(omega_, omega0_);
 445
 446     // Turbulent kinetic energy equation
 447     fvScalarMatrix kEqn
 448     (
 449         fvm::ddt(k_)
 450       + fvm::div(phi_, k_)
 451       + fvm::SuSp(-fvc::div(phi_), k_)
 452       - fvm::laplacian(DkEff(F1), k_)
 453      ==
 454         min(G, c1_*betaStar_*k_*omega_)
 455       - fvm::Sp(betaStar_*omega_, k_)
 456     );
 457
 458     kEqn.relax();
 459     solve(kEqn);
 460     bound(k_, k0_);
 461     bound(omega_, omega0_);
 462
 463     // Re-calculate viscosity
 464     // Fixed sqrt(2) error.  HJ, 10/Jun/2015
 465     nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
 466     nut_ = min(nut_, nuRatio()*nu());
 467     nut_.correctBoundaryConditions();
 468 }
 469
 470
 471 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 472
 473 } // End namespace RASModels
 474 } // End namespace incompressible
 475 } // End namespace Foam
 476
 477 // ************************************************************************* //