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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-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 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     y_(mesh_),
 200
 201     k_
 202     (
 203         IOobject
 204         (
 205             "k",
 206             runTime_.timeName(),
 207             mesh_,
 208             IOobject::NO_READ,
 209             IOobject::AUTO_WRITE
 210         ),
 211         autoCreateK("k", mesh_)
 212     ),
 213     omega_
 214     (
 215         IOobject
 216         (
 217             "omega",
 218             runTime_.timeName(),
 219             mesh_,
 220             IOobject::NO_READ,
 221             IOobject::AUTO_WRITE
 222         ),
 223         autoCreateOmega("omega", mesh_)
 224     ),
 225     nut_
 226     (
 227         IOobject
 228         (
 229             "nut",
 230             runTime_.timeName(),
 231             mesh_,
 232             IOobject::NO_READ,
 233             IOobject::AUTO_WRITE
 234         ),
 235         autoCreateNut("nut", mesh_)
 236     )
 237 {
 238     nut_ =
 239         a1_*k_
 240        /max
 241         (
 242             a1_*(omega_ + omegaSmall_),
 243             F2()*mag(symm(fvc::grad(U_)))
 244         );
 245     nut_.correctBoundaryConditions();
 246
 247     printCoeffs();
 248 }
 249
 250
 251 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 252
 253 tmp<volSymmTensorField> kOmegaSST::R() const
 254 {
 255     return tmp<volSymmTensorField>
 256     (
 257         new volSymmTensorField
 258         (
 259             IOobject
 260             (
 261                 "R",
 262                 runTime_.timeName(),
 263                 mesh_,
 264                 IOobject::NO_READ,
 265                 IOobject::NO_WRITE
 266             ),
 267             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 268             k_.boundaryField().types()
 269         )
 270     );
 271 }
 272
 273
 274 tmp<volSymmTensorField> kOmegaSST::devReff() const
 275 {
 276     return tmp<volSymmTensorField>
 277     (
 278         new volSymmTensorField
 279         (
 280             IOobject
 281             (
 282                 "devRhoReff",
 283                 runTime_.timeName(),
 284                 mesh_,
 285                 IOobject::NO_READ,
 286                 IOobject::NO_WRITE
 287             ),
 288            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 289         )
 290     );
 291 }
 292
 293
 294 tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
 295 {
 296     return
 297     (
 298       - fvm::laplacian(nuEff(), U)
 299       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 300     );
 301 }
 302
 303
 304 bool kOmegaSST::read()
 305 {
 306     if (RASModel::read())
 307     {
 308         alphaK1_.readIfPresent(coeffDict());
 309         alphaK2_.readIfPresent(coeffDict());
 310         alphaOmega1_.readIfPresent(coeffDict());
 311         alphaOmega2_.readIfPresent(coeffDict());
 312         gamma1_.readIfPresent(coeffDict());
 313         gamma2_.readIfPresent(coeffDict());
 314         beta1_.readIfPresent(coeffDict());
 315         beta2_.readIfPresent(coeffDict());
 316         betaStar_.readIfPresent(coeffDict());
 317         a1_.readIfPresent(coeffDict());
 318         c1_.readIfPresent(coeffDict());
 319
 320         return true;
 321     }
 322     else
 323     {
 324         return false;
 325     }
 326 }
 327
 328
 329 void kOmegaSST::correct()
 330 {
 331     RASModel::correct();
 332
 333     if (!turbulence_)
 334     {
 335         return;
 336     }
 337
 338     if (mesh_.changing())
 339     {
 340         y_.correct();
 341     }
 342
 343     volScalarField S2 = magSqr(symm(fvc::grad(U_)));
 344     volScalarField G("RASModel::G", nut_*2*S2);
 345
 346     // Update omega and G at the wall
 347     omega_.boundaryField().updateCoeffs();
 348
 349     volScalarField F1 = this->F1
 350     (
 351         (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
 352     );
 353
 354     // Epsilon diffusion correction
 355     surfaceScalarField CDkPhiOmega
 356     (
 357         "CDkPhiOmega",
 358         (2*alphaOmega2_)
 359        *fvc::interpolate(F1 - scalar(1))
 360        /fvc::interpolate(omega_)
 361        *fvc::snGrad(k_)*mesh_.magSf()
 362     );
 363
 364     forAll (CDkPhiOmega.boundaryField(), patchi)
 365     {
 366         if (!CDkPhiOmega.boundaryField()[patchi].coupled())
 367         {
 368             CDkPhiOmega.boundaryField()[patchi] = 0.0;
 369         }
 370     }
 371
 372     // Turbulent frequency equation
 373     tmp<fvScalarMatrix> omegaEqn
 374     (
 375         fvm::ddt(omega_)
 376       + fvm::div(phi_, omega_)
 377       - fvm::Sp(fvc::div(phi_), omega_)
 378       - fvm::laplacian(DomegaEff(F1), omega_)
 379       + fvm::div(CDkPhiOmega, omega_)
 380       - fvm::Sp(fvc::div(CDkPhiOmega), omega_)
 381      ==
 382         gamma(F1)*2*S2
 383       - fvm::Sp(beta(F1)*omega_, omega_)
 384     );
 385
 386     omegaEqn().relax();
 387
 388     omegaEqn().boundaryManipulate(omega_.boundaryField());
 389
 390     solve(omegaEqn);
 391     bound(omega_, omega0_);
 392
 393     // Turbulent kinetic energy equation
 394     tmp<fvScalarMatrix> kEqn
 395     (
 396         fvm::ddt(k_)
 397       + fvm::div(phi_, k_)
 398       - fvm::Sp(fvc::div(phi_), k_)
 399       - fvm::laplacian(DkEff(F1), k_)
 400      ==
 401         min(G, c1_*betaStar_*k_*omega_)
 402       - fvm::Sp(betaStar_*omega_, k_)
 403     );
 404
 405     kEqn().relax();
 406     solve(kEqn);
 407     bound(k_, k0_);
 408
 409
 410     // Re-calculate viscosity
 411     nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
 412     nut_.correctBoundaryConditions();
 413 }
 414
 415
 416 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 417
 418 } // End namespace RASModels
 419 } // End namespace incompressible
 420 } // End namespace Foam
 421
 422 // ************************************************************************* //