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