src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.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 "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     tmp<volScalarField> CDkOmegaPlus = max
  50     (
  51         CDkOmega,
  52         dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
  53     );
  54
  55     tmp<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 tmp<volScalarField> kOmegaSST::F2() const
  73 {
  74     tmp<volScalarField> arg2 = min
  75     (
  76         max
  77         (
  78             (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
  79             scalar(500)*nu()/(sqr(y_)*omega_)
  80         ),
  81         scalar(100)
  82     );
  83
  84     return tanh(sqr(arg2));
  85 }
  86
  87
  88 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  89
  90 kOmegaSST::kOmegaSST
  91 (
  92     const volVectorField& U,
  93     const surfaceScalarField& phi,
  94     transportModel& transport,
  95     const word& turbulenceModelName,
  96     const word& modelName
  97 )
  98 :
  99     RASModel(modelName, U, phi, transport, turbulenceModelName),
 100
 101     alphaK1_
 102     (
 103         dimensioned<scalar>::lookupOrAddToDict
 104         (
 105             "alphaK1",
 106             coeffDict_,
 107             0.85034
 108         )
 109     ),
 110     alphaK2_
 111     (
 112         dimensioned<scalar>::lookupOrAddToDict
 113         (
 114             "alphaK2",
 115             coeffDict_,
 116             1.0
 117         )
 118     ),
 119     alphaOmega1_
 120     (
 121         dimensioned<scalar>::lookupOrAddToDict
 122         (
 123             "alphaOmega1",
 124             coeffDict_,
 125             0.5
 126         )
 127     ),
 128     alphaOmega2_
 129     (
 130         dimensioned<scalar>::lookupOrAddToDict
 131         (
 132             "alphaOmega2",
 133             coeffDict_,
 134             0.85616
 135         )
 136     ),
 137     gamma1_
 138     (
 139         dimensioned<scalar>::lookupOrAddToDict
 140         (
 141             "gamma1",
 142             coeffDict_,
 143             0.5532
 144         )
 145     ),
 146     gamma2_
 147     (
 148         dimensioned<scalar>::lookupOrAddToDict
 149         (
 150             "gamma2",
 151             coeffDict_,
 152             0.4403
 153         )
 154     ),
 155     beta1_
 156     (
 157         dimensioned<scalar>::lookupOrAddToDict
 158         (
 159             "beta1",
 160             coeffDict_,
 161             0.075
 162         )
 163     ),
 164     beta2_
 165     (
 166         dimensioned<scalar>::lookupOrAddToDict
 167         (
 168             "beta2",
 169             coeffDict_,
 170             0.0828
 171         )
 172     ),
 173     betaStar_
 174     (
 175         dimensioned<scalar>::lookupOrAddToDict
 176         (
 177             "betaStar",
 178             coeffDict_,
 179             0.09
 180         )
 181     ),
 182     a1_
 183     (
 184         dimensioned<scalar>::lookupOrAddToDict
 185         (
 186             "a1",
 187             coeffDict_,
 188             0.31
 189         )
 190     ),
 191     c1_
 192     (
 193         dimensioned<scalar>::lookupOrAddToDict
 194         (
 195             "c1",
 196             coeffDict_,
 197             10.0
 198         )
 199     ),
 200
 201     y_(mesh_),
 202
 203     k_
 204     (
 205         IOobject
 206         (
 207             "k",
 208             runTime_.timeName(),
 209             mesh_,
 210             IOobject::NO_READ,
 211             IOobject::AUTO_WRITE
 212         ),
 213         autoCreateK("k", mesh_)
 214     ),
 215     omega_
 216     (
 217         IOobject
 218         (
 219             "omega",
 220             runTime_.timeName(),
 221             mesh_,
 222             IOobject::NO_READ,
 223             IOobject::AUTO_WRITE
 224         ),
 225         autoCreateOmega("omega", mesh_)
 226     ),
 227     nut_
 228     (
 229         IOobject
 230         (
 231             "nut",
 232             runTime_.timeName(),
 233             mesh_,
 234             IOobject::NO_READ,
 235             IOobject::AUTO_WRITE
 236         ),
 237         autoCreateNut("nut", mesh_)
 238     )
 239 {
 240     bound(k_, kMin_);
 241     bound(omega_, omegaMin_);
 242
 243     nut_ =
 244     (
 245         a1_*k_
 246       / max
 247         (
 248             a1_*omega_,
 249             F2()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
 250         )
 251     );
 252     nut_.correctBoundaryConditions();
 253
 254     printCoeffs();
 255 }
 256
 257
 258 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 259
 260 tmp<volSymmTensorField> kOmegaSST::R() const
 261 {
 262     return tmp<volSymmTensorField>
 263     (
 264         new volSymmTensorField
 265         (
 266             IOobject
 267             (
 268                 "R",
 269                 runTime_.timeName(),
 270                 mesh_,
 271                 IOobject::NO_READ,
 272                 IOobject::NO_WRITE
 273             ),
 274             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 275             k_.boundaryField().types()
 276         )
 277     );
 278 }
 279
 280
 281 tmp<volSymmTensorField> kOmegaSST::devReff() const
 282 {
 283     return tmp<volSymmTensorField>
 284     (
 285         new volSymmTensorField
 286         (
 287             IOobject
 288             (
 289                 "devRhoReff",
 290                 runTime_.timeName(),
 291                 mesh_,
 292                 IOobject::NO_READ,
 293                 IOobject::NO_WRITE
 294             ),
 295            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 296         )
 297     );
 298 }
 299
 300
 301 tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
 302 {
 303     return
 304     (
 305       - fvm::laplacian(nuEff(), U)
 306       - fvc::div(nuEff()*dev(T(fvc::grad(U))))
 307     );
 308 }
 309
 310
 311 bool kOmegaSST::read()
 312 {
 313     if (RASModel::read())
 314     {
 315         alphaK1_.readIfPresent(coeffDict());
 316         alphaK2_.readIfPresent(coeffDict());
 317         alphaOmega1_.readIfPresent(coeffDict());
 318         alphaOmega2_.readIfPresent(coeffDict());
 319         gamma1_.readIfPresent(coeffDict());
 320         gamma2_.readIfPresent(coeffDict());
 321         beta1_.readIfPresent(coeffDict());
 322         beta2_.readIfPresent(coeffDict());
 323         betaStar_.readIfPresent(coeffDict());
 324         a1_.readIfPresent(coeffDict());
 325         c1_.readIfPresent(coeffDict());
 326
 327         return true;
 328     }
 329     else
 330     {
 331         return false;
 332     }
 333 }
 334
 335
 336 void kOmegaSST::correct()
 337 {
 338     RASModel::correct();
 339
 340     if (!turbulence_)
 341     {
 342         return;
 343     }
 344
 345     if (mesh_.changing())
 346     {
 347         y_.correct();
 348     }
 349
 350     const volScalarField S2(2*magSqr(symm(fvc::grad(U_))));
 351     volScalarField G("RASModel::G", nut_*S2);
 352
 353     // Update omega and G at the wall
 354     omega_.boundaryField().updateCoeffs();
 355
 356     const volScalarField CDkOmega
 357     (
 358         (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
 359     );
 360
 361     const volScalarField F1(this->F1(CDkOmega));
 362
 363     // Turbulent frequency equation
 364     tmp<fvScalarMatrix> omegaEqn
 365     (
 366         fvm::ddt(omega_)
 367       + fvm::div(phi_, omega_)
 368       - fvm::Sp(fvc::div(phi_), omega_)
 369       - fvm::laplacian(DomegaEff(F1), omega_)
 370      ==
 371         gamma(F1)*S2
 372       - fvm::Sp(beta(F1)*omega_, omega_)
 373       - fvm::SuSp
 374         (
 375             (F1 - scalar(1))*CDkOmega/omega_,
 376             omega_
 377         )
 378     );
 379
 380     omegaEqn().relax();
 381
 382     omegaEqn().boundaryManipulate(omega_.boundaryField());
 383
 384     solve(omegaEqn);
 385     bound(omega_, omegaMin_);
 386
 387     // Turbulent kinetic energy equation
 388     tmp<fvScalarMatrix> kEqn
 389     (
 390         fvm::ddt(k_)
 391       + fvm::div(phi_, k_)
 392       - fvm::Sp(fvc::div(phi_), k_)
 393       - fvm::laplacian(DkEff(F1), k_)
 394      ==
 395         min(G, c1_*betaStar_*k_*omega_)
 396       - fvm::Sp(betaStar_*omega_, k_)
 397     );
 398
 399     kEqn().relax();
 400     solve(kEqn);
 401     bound(k_, kMin_);
 402
 403
 404     // Re-calculate viscosity
 405     nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
 406     nut_.correctBoundaryConditions();
 407 }
 408
 409
 410 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 411
 412 } // End namespace RASModels
 413 } // End namespace incompressible
 414 } // End namespace Foam
 415
 416 // ************************************************************************* //