src/turbulenceModels/compressible/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 "addToRunTimeSelectionTable.H"
  28
  29 #include "backwardsCompatibilityWallFunctions.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     const 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     Prt_
 137     (
 138         dimensioned<scalar>::lookupOrAddToDict
 139         (
 140             "Prt",
 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     y_(mesh_),
 210
 211     k_
 212     (
 213         IOobject
 214         (
 215             "k",
 216             runTime_.timeName(),
 217             mesh_,
 218             IOobject::NO_READ,
 219             IOobject::AUTO_WRITE
 220         ),
 221         autoCreateK("k", mesh_)
 222     ),
 223     omega_
 224     (
 225         IOobject
 226         (
 227             "omega",
 228             runTime_.timeName(),
 229             mesh_,
 230             IOobject::NO_READ,
 231             IOobject::AUTO_WRITE
 232         ),
 233         autoCreateOmega("omega", mesh_)
 234     ),
 235     mut_
 236     (
 237         IOobject
 238         (
 239             "mut",
 240             runTime_.timeName(),
 241             mesh_,
 242             IOobject::NO_READ,
 243             IOobject::AUTO_WRITE
 244         ),
 245         autoCreateMut("mut", mesh_)
 246     ),
 247     alphat_
 248     (
 249         IOobject
 250         (
 251             "alphat",
 252             runTime_.timeName(),
 253             mesh_,
 254             IOobject::NO_READ,
 255             IOobject::AUTO_WRITE
 256         ),
 257         autoCreateAlphat("alphat", mesh_)
 258     )
 259 {
 260     bound(omega_, omega0_);
 261
 262     mut_ =
 263     (
 264         a1_*rho_*k_
 265       / max
 266         (
 267             a1_*omega_,
 268             F2()*sqrt(2.0*magSqr(symm(fvc::grad(U_))))
 269         )
 270     );
 271     mut_.correctBoundaryConditions();
 272
 273     alphat_ = mut_/Prt_;
 274     alphat_.correctBoundaryConditions();
 275
 276     printCoeffs();
 277 }
 278
 279
 280 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 281
 282 tmp<volSymmTensorField> kOmegaSST::R() const
 283 {
 284     return tmp<volSymmTensorField>
 285     (
 286         new volSymmTensorField
 287         (
 288             IOobject
 289             (
 290                 "R",
 291                 runTime_.timeName(),
 292                 mesh_,
 293                 IOobject::NO_READ,
 294                 IOobject::NO_WRITE
 295             ),
 296             ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
 297             k_.boundaryField().types()
 298         )
 299     );
 300 }
 301
 302
 303 tmp<volSymmTensorField> kOmegaSST::devRhoReff() const
 304 {
 305     return tmp<volSymmTensorField>
 306     (
 307         new volSymmTensorField
 308         (
 309             IOobject
 310             (
 311                 "devRhoReff",
 312                 runTime_.timeName(),
 313                 mesh_,
 314                 IOobject::NO_READ,
 315                 IOobject::NO_WRITE
 316             ),
 317             -muEff()*dev(twoSymm(fvc::grad(U_)))
 318         )
 319     );
 320 }
 321
 322
 323 tmp<fvVectorMatrix> kOmegaSST::divDevRhoReff(volVectorField& U) const
 324 {
 325     return
 326     (
 327       - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
 328     );
 329 }
 330
 331
 332 bool kOmegaSST::read()
 333 {
 334     if (RASModel::read())
 335     {
 336         alphaK1_.readIfPresent(coeffDict());
 337         alphaK2_.readIfPresent(coeffDict());
 338         alphaOmega1_.readIfPresent(coeffDict());
 339         alphaOmega2_.readIfPresent(coeffDict());
 340         Prt_.readIfPresent(coeffDict());
 341         gamma1_.readIfPresent(coeffDict());
 342         gamma2_.readIfPresent(coeffDict());
 343         beta1_.readIfPresent(coeffDict());
 344         beta2_.readIfPresent(coeffDict());
 345         betaStar_.readIfPresent(coeffDict());
 346         a1_.readIfPresent(coeffDict());
 347         c1_.readIfPresent(coeffDict());
 348
 349         return true;
 350     }
 351     else
 352     {
 353         return false;
 354     }
 355 }
 356
 357
 358 void kOmegaSST::correct()
 359 {
 360     if (!turbulence_)
 361     {
 362         // Re-calculate viscosity
 363         mut_ =
 364             a1_*rho_*k_
 365            /max(a1_*omega_, F2()*sqrt(2.0*magSqr(symm(fvc::grad(U_)))));
 366         mut_.correctBoundaryConditions();
 367
 368         // Re-calculate thermal diffusivity
 369         alphat_ = mut_/Prt_;
 370         alphat_.correctBoundaryConditions();
 371
 372         return;
 373     }
 374
 375     RASModel::correct();
 376
 377     volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
 378
 379     if (mesh_.changing())
 380     {
 381         y_.correct();
 382     }
 383
 384     if (mesh_.moving())
 385     {
 386         divU += fvc::div(mesh_.phi());
 387     }
 388
 389     tmp<volTensorField> tgradU = fvc::grad(U_);
 390     volScalarField S2 = magSqr(symm(tgradU()));
 391     volScalarField GbyMu = (tgradU() && dev(twoSymm(tgradU())));
 392     volScalarField G("RASModel::G", mut_*GbyMu);
 393     tgradU.clear();
 394
 395     // Update omega and G at the wall
 396     omega_.boundaryField().updateCoeffs();
 397
 398     volScalarField CDkOmega =
 399         (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
 400
 401     volScalarField F1 = this->F1(CDkOmega);
 402     volScalarField rhoGammaF1 = rho_*gamma(F1);
 403
 404     // Turbulent frequency equation
 405     tmp<fvScalarMatrix> omegaEqn
 406     (
 407         fvm::ddt(rho_, omega_)
 408       + fvm::div(phi_, omega_)
 409       - fvm::laplacian(DomegaEff(F1), omega_)
 410      ==
 411         rhoGammaF1*GbyMu
 412       - fvm::SuSp((2.0/3.0)*rhoGammaF1*divU, omega_)
 413       - fvm::Sp(rho_*beta(F1)*omega_, omega_)
 414       - fvm::SuSp
 415         (
 416             rho_*(F1 - scalar(1))*CDkOmega/omega_,
 417             omega_
 418         )
 419     );
 420
 421     omegaEqn().relax();
 422
 423     omegaEqn().boundaryManipulate(omega_.boundaryField());
 424
 425     solve(omegaEqn);
 426     bound(omega_, omega0_);
 427
 428     // Turbulent kinetic energy equation
 429     tmp<fvScalarMatrix> kEqn
 430     (
 431         fvm::ddt(rho_, k_)
 432       + fvm::div(phi_, k_)
 433       - fvm::laplacian(DkEff(F1), k_)
 434      ==
 435         min(G, (c1_*betaStar_)*rho_*k_*omega_)
 436       - fvm::SuSp(2.0/3.0*rho_*divU, k_)
 437       - fvm::Sp(rho_*betaStar_*omega_, k_)
 438     );
 439
 440     kEqn().relax();
 441     solve(kEqn);
 442     bound(k_, k0_);
 443
 444
 445     // Re-calculate viscosity
 446     mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(2.0*S2));
 447     mut_.correctBoundaryConditions();
 448
 449     // Re-calculate thermal diffusivity
 450     alphat_ = mut_/Prt_;
 451     alphat_.correctBoundaryConditions();
 452 }
 453
 454
 455 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 456
 457 } // End namespace RASModels
 458 } // End namespace compressible
 459 } // End namespace Foam
 460
 461 // ************************************************************************* //