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