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