src/turbulenceModels/incompressible/RAS/coupledKOmegaSST/coupledKOmegaSST.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "coupledKOmegaSST.H"
  27 #include "fvBlockMatrix.H"
  28 #include "addToRunTimeSelectionTable.H"
  29
  30 #include "backwardsCompatibilityWallFunctions.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36 namespace incompressible
  37 {
  38 namespace RASModels
  39 {
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 defineTypeNameAndDebug(coupledKOmegaSST, 0);
  44 addToRunTimeSelectionTable(RASModel, coupledKOmegaSST, dictionary);
  45
  46 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  47
  48 tmp<volScalarField> coupledKOmegaSST::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)*nu()/(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> coupledKOmegaSST::F2() const
  74 {
  75     volScalarField arg2 = min
  76     (
  77         max
  78         (
  79             (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
  80             scalar(500)*nu()/(sqr(y_)*omega_)
  81         ),
  82         scalar(100)
  83     );
  84
  85     return tanh(sqr(arg2));
  86 }
  87
  88
  89 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  90
  91 coupledKOmegaSST::coupledKOmegaSST
  92 (
  93     const volVectorField& U,
  94     const surfaceScalarField& phi,
  95     transportModel& lamTransportModel
  96 )
  97 :
  98     RASModel(typeName, U, phi, lamTransportModel),
  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     gamma1_
 137     (
 138         dimensioned<scalar>::lookupOrAddToDict
 139         (
 140             "gamma1",
 141             coeffDict_,
 142             0.5532
 143         )
 144     ),
 145     gamma2_
 146     (
 147         dimensioned<scalar>::lookupOrAddToDict
 148         (
 149             "gamma2",
 150             coeffDict_,
 151             0.4403
 152         )
 153     ),
 154     beta1_
 155     (
 156         dimensioned<scalar>::lookupOrAddToDict
 157         (
 158             "beta1",
 159             coeffDict_,
 160             0.075
 161         )
 162     ),
 163     beta2_
 164     (
 165         dimensioned<scalar>::lookupOrAddToDict
 166         (
 167             "beta2",
 168             coeffDict_,
 169             0.0828
 170         )
 171     ),
 172     betaStar_
 173     (
 174         dimensioned<scalar>::lookupOrAddToDict
 175         (
 176             "betaStar",
 177             coeffDict_,
 178             0.09
 179         )
 180     ),
 181     a1_
 182     (
 183         dimensioned<scalar>::lookupOrAddToDict
 184         (
 185             "a1",
 186             coeffDict_,
 187             0.31
 188         )
 189     ),
 190     c1_
 191     (
 192         dimensioned<scalar>::lookupOrAddToDict
 193         (
 194             "c1",
 195             coeffDict_,
 196             10.0
 197         )
 198     ),
 199
 200     y_(mesh_),
 201
 202     k_
 203     (
 204         IOobject
 205         (
 206             "k",
 207             runTime_.timeName(),
 208             U_.db(),
 209             IOobject::NO_READ,
 210             IOobject::AUTO_WRITE
 211         ),
 212         autoCreateK("k", mesh_, U_.db())
 213     ),
 214     omega_
 215     (
 216         IOobject
 217         (
 218             "omega",
 219             runTime_.timeName(),
 220             U_.db(),
 221             IOobject::NO_READ,
 222             IOobject::AUTO_WRITE
 223         ),
 224         autoCreateOmega("omega", mesh_, U_.db())
 225     ),
 226     nut_
 227     (
 228         IOobject
 229         (
 230             "nut",
 231             runTime_.timeName(),
 232             U_.db(),
 233             IOobject::NO_READ,
 234             IOobject::AUTO_WRITE
 235         ),
 236         autoCreateNut("nut", mesh_, U_.db())
 237     ),
 238     kOmega_
 239     (
 240         IOobject
 241         (
 242             "kOmega",
 243             runTime_.timeName(),
 244             U_.db(),
 245             IOobject::NO_READ,
 246             IOobject::NO_WRITE
 247         ),
 248         mesh_,
 249         dimensionedVector2("zero", dimless, vector2::zero)
 250     )
 251 {
 252     bound(omega_, omega0_);
 253
 254     nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2.0)*mag(symm(fvc::grad(U_))));
 255     nut_.correctBoundaryConditions();
 256
 257     printCoeffs();
 258 }
 259
 260
 261 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 262
 263 tmp<volSymmTensorField> coupledKOmegaSST::R() const
 264 {
 265     return tmp<volSymmTensorField>
 266     (
 267         new volSymmTensorField
 268         (
 269             IOobject
 270             (
 271                 "R",
 272                 runTime_.timeName(),
 273                 U_.db(),
 274                 IOobject::NO_READ,
 275                 IOobject::NO_WRITE
 276             ),
 277             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 278             k_.boundaryField().types()
 279         )
 280     );
 281 }
 282
 283
 284 tmp<volSymmTensorField> coupledKOmegaSST::devReff() const
 285 {
 286     return tmp<volSymmTensorField>
 287     (
 288         new volSymmTensorField
 289         (
 290             IOobject
 291             (
 292                 "devRhoReff",
 293                 runTime_.timeName(),
 294                 U_.db(),
 295                 IOobject::NO_READ,
 296                 IOobject::NO_WRITE
 297             ),
 298            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 299         )
 300     );
 301 }
 302
 303
 304 tmp<fvVectorMatrix> coupledKOmegaSST::divDevReff(volVectorField& U) const
 305 {
 306     return
 307     (
 308       - fvm::laplacian(nuEff(), U)
 309       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 310     );
 311 }
 312
 313
 314 bool coupledKOmegaSST::read()
 315 {
 316     if (RASModel::read())
 317     {
 318         alphaK1_.readIfPresent(coeffDict());
 319         alphaK2_.readIfPresent(coeffDict());
 320         alphaOmega1_.readIfPresent(coeffDict());
 321         alphaOmega2_.readIfPresent(coeffDict());
 322         gamma1_.readIfPresent(coeffDict());
 323         gamma2_.readIfPresent(coeffDict());
 324         beta1_.readIfPresent(coeffDict());
 325         beta2_.readIfPresent(coeffDict());
 326         betaStar_.readIfPresent(coeffDict());
 327         a1_.readIfPresent(coeffDict());
 328         c1_.readIfPresent(coeffDict());
 329
 330         return true;
 331     }
 332     else
 333     {
 334         return false;
 335     }
 336 }
 337
 338
 339 void coupledKOmegaSST::correct()
 340 {
 341     // Bound in case of topological change
 342     // HJ, 22/Aug/2007
 343     if (mesh_.changing())
 344     {
 345         bound(k_, k0_);
 346         bound(omega_, omega0_);
 347     }
 348
 349     RASModel::correct();
 350
 351     if (!turbulence_)
 352     {
 353         return;
 354     }
 355
 356     if (mesh_.changing())
 357     {
 358         y_.correct();
 359     }
 360
 361     volScalarField S2 = magSqr(symm(fvc::grad(U_)));
 362     volScalarField G("RASModel::G", nut_*2*S2);
 363
 364     // Make coupled matrix
 365     fvBlockMatrix<vector2> kOmegaEqn(kOmega_);
 366
 367     // Update omega and G at the wall
 368     omega_.boundaryField().updateCoeffs();
 369
 370     volScalarField CDkOmega =
 371         (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
 372
 373     volScalarField F1 = this->F1(CDkOmega);
 374
 375     Info<< "Coupled k-omega" << endl;
 376
 377     // Turbulent frequency equation
 378     {
 379         fvScalarMatrix omegaEqn
 380         (
 381             fvm::ddt(omega_)
 382           + fvm::div(phi_, omega_)
 383           + fvm::SuSp(-fvc::div(phi_), omega_)
 384           - fvm::laplacian(DomegaEff(F1), omega_)
 385           + fvm::SuSp
 386             (
 387                 beta(F1)*omega_
 388               + (F1 - scalar(1))*CDkOmega/omega_,
 389                 omega_
 390             )
 391          ==
 392             gamma(F1)*2*S2
 393         );
 394
 395         omegaEqn.relax();
 396         omegaEqn.completeAssembly();
 397
 398         kOmegaEqn.insertEquation(1, omegaEqn);
 399
 400         // Add coupling term
 401         volScalarField coupling
 402         (
 403             "coupling",
 404             -gamma(F1)*2*S2/k_
 405         );
 406         scalarField& couplingIn = coupling.internalField();
 407
 408         // Eliminate coupling in wall function cells
 409         labelList eliminated = omegaEqn.eliminatedEqns().toc();
 410
 411         forAll (eliminated, cellI)
 412         {
 413             couplingIn[eliminated[cellI]] *= 0;
 414         }
 415
 416         // Insert coupling
 417         kOmegaEqn.insertEquationCoupling(1, 0, coupling);
 418     }
 419
 420     // Turbulent kinetic energy equation
 421     {
 422         fvScalarMatrix kEqn
 423         (
 424             fvm::ddt(k_)
 425           + fvm::div(phi_, k_)
 426           + fvm::SuSp(-fvc::div(phi_), k_)
 427           - fvm::laplacian(DkEff(F1), k_)
 428           + fvm::SuSp
 429             (
 430                 betaStar_*omega_
 431               - min(G/k_, c1_*betaStar_*omega_),
 432                 k_
 433             )
 434         );
 435
 436         kEqn.relax();
 437
 438         kOmegaEqn.insertEquation(0, kEqn);
 439     }
 440
 441     // Update source coupling: coupling terms eliminated from source
 442     kOmegaEqn.updateSourceCoupling();
 443
 444     kOmegaEqn.solve();
 445
 446     // Retrieve solution
 447     kOmegaEqn.retrieveSolution(0, k_.internalField());
 448     kOmegaEqn.retrieveSolution(1, omega_.internalField());
 449
 450     k_.correctBoundaryConditions();
 451     omega_.correctBoundaryConditions();
 452
 453     bound(k_, k0_);
 454     bound(omega_, omega0_);
 455
 456     // Re-calculate viscosity
 457     nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2*S2));
 458     nut_ = min(nut_, nuRatio()*nu());
 459     nut_.correctBoundaryConditions();
 460 }
 461
 462
 463 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 464
 465 } // End namespace RASModels
 466 } // End namespace incompressible
 467 } // End namespace Foam
 468
 469 // ************************************************************************* //