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