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