src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "kOmegaSST.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "wallFvPatch.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     omega0_("omega0", dimless/dimTime, SMALL),
 200     omegaSmall_("omegaSmall", dimless/dimTime, SMALL),
 201
 202     Cmu_
 203     (
 204         dimensioned<scalar>::lookupOrAddToDict
 205         (
 206             "Cmu",
 207             coeffDict_,
 208             0.09
 209         )
 210     ),
 211
 212     y_(mesh_),
 213
 214     k_
 215     (
 216         IOobject
 217         (
 218             "k",
 219             runTime_.timeName(),
 220             mesh_,
 221             IOobject::MUST_READ,
 222             IOobject::AUTO_WRITE
 223         ),
 224         mesh_
 225     ),
 226
 227     omega_
 228     (
 229         IOobject
 230         (
 231             "omega",
 232             runTime_.timeName(),
 233             mesh_,
 234             IOobject::MUST_READ,
 235             IOobject::AUTO_WRITE
 236         ),
 237         mesh_
 238     ),
 239
 240     nut_(a1_*k_/max(a1_*(omega_ + omegaSmall_), F2()*mag(symm(fvc::grad(U_)))))
 241 {
 242 #   include "kOmegaWallViscosityI.H"
 243
 244     printCoeffs();
 245 }
 246
 247
 248 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 249
 250 tmp<volSymmTensorField> kOmegaSST::R() const
 251 {
 252     return tmp<volSymmTensorField>
 253     (
 254         new volSymmTensorField
 255         (
 256             IOobject
 257             (
 258                 "R",
 259                 runTime_.timeName(),
 260                 mesh_,
 261                 IOobject::NO_READ,
 262                 IOobject::NO_WRITE
 263             ),
 264             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 265             k_.boundaryField().types()
 266         )
 267     );
 268 }
 269
 270
 271 tmp<volSymmTensorField> kOmegaSST::devReff() const
 272 {
 273     return tmp<volSymmTensorField>
 274     (
 275         new volSymmTensorField
 276         (
 277             IOobject
 278             (
 279                 "devRhoReff",
 280                 runTime_.timeName(),
 281                 mesh_,
 282                 IOobject::NO_READ,
 283                 IOobject::NO_WRITE
 284             ),
 285            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 286         )
 287     );
 288 }
 289
 290
 291 tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
 292 {
 293     return
 294     (
 295       - fvm::laplacian(nuEff(), U)
 296       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 297     );
 298 }
 299
 300
 301 bool kOmegaSST::read()
 302 {
 303     if (RASModel::read())
 304     {
 305         alphaK1_.readIfPresent(coeffDict_);
 306         alphaK2_.readIfPresent(coeffDict_);
 307         alphaOmega1_.readIfPresent(coeffDict_);
 308         alphaOmega2_.readIfPresent(coeffDict_);
 309         gamma1_.readIfPresent(coeffDict_);
 310         gamma2_.readIfPresent(coeffDict_);
 311         beta1_.readIfPresent(coeffDict_);
 312         beta2_.readIfPresent(coeffDict_);
 313         betaStar_.readIfPresent(coeffDict_);
 314         a1_.readIfPresent(coeffDict_);
 315         c1_.readIfPresent(coeffDict_);
 316         Cmu_.readIfPresent(coeffDict_);
 317
 318         return true;
 319     }
 320     else
 321     {
 322         return false;
 323     }
 324 }
 325
 326
 327 void kOmegaSST::correct()
 328 {
 329     transportModel_.correct();
 330
 331     if (!turbulence_)
 332     {
 333         return;
 334     }
 335
 336     RASModel::correct();
 337
 338     if (mesh_.changing())
 339     {
 340         y_.correct();
 341     }
 342
 343     volScalarField S2 = magSqr(symm(fvc::grad(U_)));
 344     volScalarField G = nut_*2*S2;
 345
 346 #   include "kOmegaWallFunctionsI.H"
 347
 348     volScalarField CDkOmega =
 349         (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
 350
 351     volScalarField F1 = this->F1(CDkOmega);
 352
 353     // Turbulent frequency equation
 354     tmp<fvScalarMatrix> omegaEqn
 355     (
 356         fvm::ddt(omega_)
 357       + fvm::div(phi_, omega_)
 358       - fvm::Sp(fvc::div(phi_), omega_)
 359       - fvm::laplacian(DomegaEff(F1), omega_)
 360      ==
 361         gamma(F1)*2*S2
 362       - fvm::Sp(beta(F1)*omega_, omega_)
 363       - fvm::SuSp
 364         (
 365             (F1 - scalar(1))*CDkOmega/omega_,
 366             omega_
 367         )
 368     );
 369
 370     omegaEqn().relax();
 371
 372 #   include "wallOmegaI.H"
 373
 374     solve(omegaEqn);
 375     bound(omega_, omega0_);
 376
 377     // Turbulent kinetic energy equation
 378     tmp<fvScalarMatrix> kEqn
 379     (
 380         fvm::ddt(k_)
 381       + fvm::div(phi_, k_)
 382       - fvm::Sp(fvc::div(phi_), k_)
 383       - fvm::laplacian(DkEff(F1), k_)
 384      ==
 385         min(G, c1_*betaStar_*k_*omega_)
 386       - fvm::Sp(betaStar_*omega_, k_)
 387     );
 388
 389     kEqn().relax();
 390     solve(kEqn);
 391     bound(k_, k0_);
 392
 393
 394     // Re-calculate viscosity
 395     nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
 396
 397 #   include "kOmegaWallViscosityI.H"
 398
 399 }
 400
 401
 402 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 403
 404 } // End namespace RASModels
 405 } // End namespace incompressible
 406 } // End namespace Foam
 407
 408 // ************************************************************************* //