src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.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 "LienCubicKELowRe.H"
  27 #include "wallFvPatch.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(LienCubicKELowRe, 0);
  44 addToRunTimeSelectionTable(RASModel, LienCubicKELowRe, dictionary);
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 LienCubicKELowRe::LienCubicKELowRe
  49 (
  50     const volVectorField& U,
  51     const surfaceScalarField& phi,
  52     transportModel& lamTransportModel
  53 )
  54 :
  55     RASModel(typeName, U, phi, lamTransportModel),
  56
  57     C1_
  58     (
  59         dimensionedScalar::lookupOrAddToDict
  60         (
  61             "C1",
  62             coeffDict_,
  63             1.44
  64         )
  65     ),
  66     C2_
  67     (
  68         dimensionedScalar::lookupOrAddToDict
  69         (
  70             "C2",
  71             coeffDict_,
  72             1.92
  73         )
  74     ),
  75     sigmak_
  76     (
  77         dimensionedScalar::lookupOrAddToDict
  78         (
  79             "sigmak",
  80             coeffDict_,
  81             1.0
  82         )
  83     ),
  84     sigmaEps_
  85     (
  86         dimensionedScalar::lookupOrAddToDict
  87         (
  88             "sigmaEps",
  89             coeffDict_,
  90             1.3
  91         )
  92     ),
  93     A1_
  94     (
  95         dimensionedScalar::lookupOrAddToDict
  96         (
  97             "A1",
  98             coeffDict_,
  99             1.25
 100         )
 101     ),
 102     A2_
 103     (
 104         dimensionedScalar::lookupOrAddToDict
 105         (
 106             "A2",
 107             coeffDict_,
 108             1000.0
 109         )
 110     ),
 111     Ctau1_
 112     (
 113         dimensionedScalar::lookupOrAddToDict
 114         (
 115             "Ctau1",
 116             coeffDict_,
 117             -4.0
 118         )
 119     ),
 120     Ctau2_
 121     (
 122         dimensionedScalar::lookupOrAddToDict
 123         (
 124             "Ctau2",
 125             coeffDict_,
 126             13.0
 127         )
 128     ),
 129     Ctau3_
 130     (
 131         dimensionedScalar::lookupOrAddToDict
 132         (
 133             "Ctau3",
 134             coeffDict_,
 135             -2.0
 136         )
 137     ),
 138     alphaKsi_
 139     (
 140         dimensionedScalar::lookupOrAddToDict
 141         (
 142             "alphaKsi",
 143             coeffDict_,
 144             0.9
 145         )
 146     ),
 147     CmuWall_
 148     (
 149         dimensionedScalar::lookupOrAddToDict
 150         (
 151             "Cmu",
 152             coeffDict_,
 153             0.09
 154         )
 155     ),
 156     kappa_
 157     (
 158         dimensionedScalar::lookupOrAddToDict
 159         (
 160             "kappa",
 161             coeffDict_,
 162             0.41
 163         )
 164     ),
 165     Am_
 166     (
 167         dimensionedScalar::lookupOrAddToDict
 168         (
 169             "Am",
 170             coeffDict_,
 171             0.016
 172         )
 173     ),
 174     Aepsilon_
 175     (
 176         dimensionedScalar::lookupOrAddToDict
 177         (
 178             "Aepsilon",
 179             coeffDict_,
 180             0.263
 181         )
 182     ),
 183     Amu_
 184     (
 185         dimensionedScalar::lookupOrAddToDict
 186         (
 187             "Amu",
 188             coeffDict_,
 189             0.00222
 190         )
 191     ),
 192
 193     k_
 194     (
 195         IOobject
 196         (
 197             "k",
 198             runTime_.timeName(),
 199             U_.db(),
 200             IOobject::MUST_READ,
 201             IOobject::AUTO_WRITE
 202         ),
 203         mesh_
 204     ),
 205
 206     epsilon_
 207     (
 208         IOobject
 209         (
 210             "epsilon",
 211             runTime_.timeName(),
 212             U_.db(),
 213             IOobject::MUST_READ,
 214             IOobject::AUTO_WRITE
 215         ),
 216         mesh_
 217     ),
 218
 219     y_(mesh_),
 220
 221     gradU_(fvc::grad(U)),
 222     eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
 223     ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
 224     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
 225     fEta_(A2_ + pow(eta_, 3.0)),
 226
 227     C5viscosity_
 228     (
 229         -2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
 230        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()))
 231     ),
 232
 233     yStar_(sqrt(k_)*y_/nu() + SMALL),
 234
 235     nut_
 236     (
 237         IOobject
 238         (
 239             "nut",
 240             runTime_.timeName(),
 241             U_.db(),
 242             IOobject::NO_READ,
 243             IOobject::AUTO_WRITE
 244         ),
 245         autoCreateLowReNut("nut", mesh_)
 246     ),
 247
 248     nonlinearStress_
 249     (
 250         "nonlinearStress",
 251         symm
 252         (
 253             // quadratic terms
 254             pow(k_, 3.0)/sqr(epsilon_)
 255            *(
 256                 Ctau1_/fEta_
 257                *(
 258                     (gradU_ & gradU_)
 259                   + (gradU_ & gradU_)().T()
 260                 )
 261               + Ctau2_/fEta_*(gradU_ & gradU_.T())
 262               + Ctau3_/fEta_*(gradU_.T() & gradU_)
 263             )
 264             // cubic term C4
 265           - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
 266            *pow(Cmu_, 3.0)
 267            *(
 268                 ((gradU_ & gradU_) & gradU_.T())
 269               + ((gradU_ & gradU_.T()) & gradU_.T())
 270               - ((gradU_.T() & gradU_) & gradU_)
 271               - ((gradU_.T() & gradU_.T()) & gradU_)
 272             )
 273             // cubic term C5, explicit part
 274           + min
 275             (
 276                 C5viscosity_,
 277                 dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
 278             )*gradU_
 279         )
 280     )
 281 {
 282     nut_ = Cmu_
 283        *(
 284             scalar(1) - exp(-Am_*yStar_))
 285            /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL
 286         )
 287        *sqr(k_)/(epsilon_ + epsilonSmall_)
 288         // cubic term C5, implicit part
 289       + max
 290         (
 291             C5viscosity_,
 292             dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
 293         );
 294
 295     nut_ = min(nut_, nuRatio()*nu());
 296     nut_.correctBoundaryConditions();
 297
 298     printCoeffs();
 299 }
 300
 301
 302 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 303
 304 tmp<volSymmTensorField> LienCubicKELowRe::R() const
 305 {
 306     return tmp<volSymmTensorField>
 307     (
 308         new volSymmTensorField
 309         (
 310             IOobject
 311             (
 312                 "R",
 313                 runTime_.timeName(),
 314                 U_.db(),
 315                 IOobject::NO_READ,
 316                 IOobject::NO_WRITE
 317             ),
 318             ((2.0/3.0)*I)*k_ - nut_*twoSymm(gradU_) + nonlinearStress_,
 319             k_.boundaryField().types()
 320         )
 321     );
 322 }
 323
 324
 325 tmp<volSymmTensorField> LienCubicKELowRe::devReff() const
 326 {
 327     return tmp<volSymmTensorField>
 328     (
 329         new volSymmTensorField
 330         (
 331             IOobject
 332             (
 333                 "devRhoReff",
 334                 runTime_.timeName(),
 335                 U_.db(),
 336                 IOobject::NO_READ,
 337                 IOobject::NO_WRITE
 338             ),
 339            -nuEff()*dev(twoSymm(fvc::grad(U_))) + nonlinearStress_
 340         )
 341     );
 342 }
 343
 344
 345 tmp<fvVectorMatrix> LienCubicKELowRe::divDevReff(volVectorField& U) const
 346 {
 347     return
 348     (
 349         fvc::div(nonlinearStress_)
 350       - fvm::laplacian(nuEff(), U)
 351       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 352     );
 353 }
 354
 355
 356 bool LienCubicKELowRe::read()
 357 {
 358     if (RASModel::read())
 359     {
 360         C1_.readIfPresent(coeffDict());
 361         C2_.readIfPresent(coeffDict());
 362         sigmak_.readIfPresent(coeffDict());
 363         sigmaEps_.readIfPresent(coeffDict());
 364         A1_.readIfPresent(coeffDict());
 365         A2_.readIfPresent(coeffDict());
 366         Ctau1_.readIfPresent(coeffDict());
 367         Ctau2_.readIfPresent(coeffDict());
 368         Ctau3_.readIfPresent(coeffDict());
 369         alphaKsi_.readIfPresent(coeffDict());
 370         CmuWall_.readIfPresent(coeffDict());
 371         kappa_.readIfPresent(coeffDict());
 372         Am_.readIfPresent(coeffDict());
 373         Aepsilon_.readIfPresent(coeffDict());
 374         Amu_.readIfPresent(coeffDict());
 375
 376         return true;
 377     }
 378     else
 379     {
 380         return false;
 381     }
 382 }
 383
 384
 385 void LienCubicKELowRe::correct()
 386 {
 387     // Bound in case of topological change
 388     // HJ, 22/Aug/2007
 389     if (mesh_.changing())
 390     {
 391         bound(k_, k0_);
 392         bound(epsilon_, epsilon0_);
 393     }
 394
 395     RASModel::correct();
 396
 397     if (!turbulence_)
 398     {
 399         return;
 400     }
 401
 402     if (mesh_.changing())
 403     {
 404         y_.correct();
 405     }
 406
 407     gradU_ = fvc::grad(U_);
 408
 409     // generation term
 410     volScalarField S2 = symm(gradU_) && gradU_;
 411
 412     yStar_ = sqrt(k_)*y_/nu() + SMALL;
 413     volScalarField Rt = sqr(k_)/(nu()*epsilon_);
 414
 415     volScalarField fMu =
 416         (scalar(1) - exp(-Am_*yStar_))
 417        /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL);
 418
 419     volScalarField f2 = scalar(1) - 0.3*exp(-sqr(Rt));
 420
 421     volScalarField G =
 422         Cmu_*fMu*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_);
 423
 424     // Dissipation equation
 425     tmp<fvScalarMatrix> epsEqn
 426     (
 427         fvm::ddt(epsilon_)
 428       + fvm::div(phi_, epsilon_)
 429       + fvm::SuSp(-fvc::div(phi_), epsilon_)
 430       - fvm::laplacian(DepsilonEff(), epsilon_)
 431       ==
 432         C1_*G*epsilon_/k_
 433         // E-term
 434       + C2_*f2*pow(Cmu_, 0.75)*pow(k_, scalar(0.5))
 435        /(kappa_*y_*(scalar(1) - exp(-Aepsilon_*yStar_)))
 436        *exp(-Amu_*sqr(yStar_))*epsilon_
 437       - fvm::Sp(C2_*f2*epsilon_/k_, epsilon_)
 438     );
 439
 440     epsEqn().relax();
 441
 442 #   include "LienCubicKELowReSetWallDissipation.H"
 443 #   include "wallDissipationI.H"
 444
 445     solve(epsEqn);
 446     bound(epsilon_, epsilon0_);
 447
 448
 449     // Turbulent kinetic energy equation
 450
 451     tmp<fvScalarMatrix> kEqn
 452     (
 453         fvm::ddt(k_)
 454       + fvm::div(phi_, k_)
 455       + fvm::SuSp(-fvc::div(phi_), k_)
 456       - fvm::laplacian(DkEff(), k_)
 457       ==
 458         G
 459       - fvm::Sp(epsilon_/k_, k_)
 460     );
 461
 462     kEqn().relax();
 463     solve(kEqn);
 464     bound(k_, k0_);
 465
 466
 467     // Re-calculate viscosity
 468
 469     eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())));
 470     ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())));
 471     Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
 472     fEta_ = A2_ + pow(eta_, 3.0);
 473
 474     C5viscosity_ =
 475       - 2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
 476        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()));
 477
 478     nut_ =
 479         Cmu_*fMu*sqr(k_)/epsilon_
 480         // C5 term, implicit
 481       + max
 482         (
 483             C5viscosity_,
 484             dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
 485         );
 486     nut_ = min(nut_, nuRatio()*nu());
 487     nut_.correctBoundaryConditions();
 488
 489     nonlinearStress_ = symm
 490     (
 491         // quadratic terms
 492         pow(k_, 3.0)/sqr(epsilon_)
 493        *(
 494             Ctau1_/fEta_
 495            *(
 496                 (gradU_ & gradU_)
 497               + (gradU_ & gradU_)().T()
 498             )
 499           + Ctau2_/fEta_*(gradU_ & gradU_.T())
 500           + Ctau3_/fEta_*(gradU_.T() & gradU_)
 501         )
 502         // cubic term C4
 503       - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
 504        *pow(Cmu_, 3.0)
 505        *(
 506             ((gradU_ & gradU_) & gradU_.T())
 507           + ((gradU_ & gradU_.T()) & gradU_.T())
 508           - ((gradU_.T() & gradU_) & gradU_)
 509           - ((gradU_.T() & gradU_.T()) & gradU_)
 510         )
 511         // cubic term C5, explicit part
 512       + min
 513         (
 514             C5viscosity_,
 515             dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
 516         )*gradU_
 517     );
 518 }
 519
 520
 521 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 522
 523 } // End namespace RASModels
 524 } // End namespace incompressible
 525 } // End namespace Foam
 526
 527 // ************************************************************************* //