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