src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.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 "LienCubicKE.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(LienCubicKE, 0);
  44 addToRunTimeSelectionTable(RASModel, LienCubicKE, dictionary);
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 LienCubicKE::LienCubicKE
  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         dimensioned<scalar>::lookupOrAddToDict
  60         (
  61             "C1",
  62             coeffDict_,
  63             1.44
  64         )
  65     ),
  66     C2_
  67     (
  68         dimensioned<scalar>::lookupOrAddToDict
  69         (
  70             "C2",
  71             coeffDict_,
  72             1.92
  73         )
  74     ),
  75     sigmak_
  76     (
  77         dimensioned<scalar>::lookupOrAddToDict
  78         (
  79             "sigmak",
  80             coeffDict_,
  81             1.0
  82         )
  83     ),
  84     sigmaEps_
  85     (
  86         dimensioned<scalar>::lookupOrAddToDict
  87         (
  88             "sigmaEps",
  89             coeffDict_,
  90             1.3
  91         )
  92     ),
  93     A1_
  94     (
  95         dimensioned<scalar>::lookupOrAddToDict
  96         (
  97             "A1",
  98             coeffDict_,
  99             1.25
 100         )
 101     ),
 102     A2_
 103     (
 104         dimensioned<scalar>::lookupOrAddToDict
 105         (
 106             "A2",
 107             coeffDict_,
 108             1000.0
 109         )
 110     ),
 111     Ctau1_
 112     (
 113         dimensioned<scalar>::lookupOrAddToDict
 114         (
 115             "Ctau1",
 116             coeffDict_,
 117             -4.0
 118         )
 119     ),
 120     Ctau2_
 121     (
 122         dimensioned<scalar>::lookupOrAddToDict
 123         (
 124             "Ctau2",
 125             coeffDict_,
 126             13.0
 127         )
 128     ),
 129     Ctau3_
 130     (
 131         dimensioned<scalar>::lookupOrAddToDict
 132         (
 133             "Ctau3",
 134             coeffDict_,
 135             -2.0
 136         )
 137     ),
 138     alphaKsi_
 139     (
 140         dimensioned<scalar>::lookupOrAddToDict
 141         (
 142             "alphaKsi",
 143             coeffDict_,
 144             0.9
 145         )
 146     ),
 147
 148     k_
 149     (
 150         IOobject
 151         (
 152             "k",
 153             runTime_.timeName(),
 154             mesh_,
 155             IOobject::NO_READ,
 156             IOobject::AUTO_WRITE
 157         ),
 158         autoCreateK("k", mesh_)
 159     ),
 160     epsilon_
 161     (
 162         IOobject
 163         (
 164             "epsilon",
 165             runTime_.timeName(),
 166             mesh_,
 167             IOobject::NO_READ,
 168             IOobject::AUTO_WRITE
 169         ),
 170         autoCreateEpsilon("epsilon", mesh_)
 171     ),
 172
 173     gradU_(fvc::grad(U)),
 174     eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
 175     ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
 176     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
 177     fEta_(A2_ + pow(eta_, 3.0)),
 178
 179     C5viscosity_
 180     (
 181       - 2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
 182        *(
 183             magSqr(gradU_ + gradU_.T())
 184           - magSqr(gradU_ - gradU_.T())
 185         )
 186     ),
 187
 188     nut_
 189     (
 190         IOobject
 191         (
 192             "nut",
 193             runTime_.timeName(),
 194             mesh_,
 195             IOobject::NO_READ,
 196             IOobject::AUTO_WRITE
 197         ),
 198         autoCreateNut("nut", mesh_)
 199     ),
 200
 201     nonlinearStress_
 202     (
 203         "nonlinearStress",
 204         // quadratic terms
 205         symm
 206         (
 207             pow(k_, 3.0)/sqr(epsilon_)
 208            *(
 209                 Ctau1_/fEta_
 210                *(
 211                     (gradU_ & gradU_)
 212                   + (gradU_ & gradU_)().T()
 213                 )
 214               + Ctau2_/fEta_*(gradU_ & gradU_.T())
 215               + Ctau3_/fEta_*(gradU_.T() & gradU_)
 216             )
 217             // cubic term C4
 218           - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
 219            *pow(Cmu_, 3.0)
 220            *(
 221                 ((gradU_ & gradU_) & gradU_.T())
 222               + ((gradU_ & gradU_.T()) & gradU_.T())
 223               - ((gradU_.T() & gradU_) & gradU_)
 224               - ((gradU_.T() & gradU_.T()) & gradU_)
 225             )
 226         )
 227     )
 228 {
 229     nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_;
 230     nut_.correctBoundaryConditions();
 231
 232     printCoeffs();
 233 }
 234
 235
 236 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 237
 238 tmp<volSymmTensorField> LienCubicKE::R() const
 239 {
 240     return tmp<volSymmTensorField>
 241     (
 242         new volSymmTensorField
 243         (
 244             IOobject
 245             (
 246                 "R",
 247                 runTime_.timeName(),
 248                 mesh_,
 249                 IOobject::NO_READ,
 250                 IOobject::NO_WRITE
 251             ),
 252             ((2.0/3.0)*I)*k_ - nut_*twoSymm(gradU_) + nonlinearStress_,
 253             k_.boundaryField().types()
 254         )
 255     );
 256 }
 257
 258
 259 tmp<volSymmTensorField> LienCubicKE::devReff() const
 260 {
 261     return tmp<volSymmTensorField>
 262     (
 263         new volSymmTensorField
 264         (
 265             IOobject
 266             (
 267                 "devRhoReff",
 268                 runTime_.timeName(),
 269                 mesh_,
 270                 IOobject::NO_READ,
 271                 IOobject::NO_WRITE
 272             ),
 273            -nuEff()*dev(twoSymm(fvc::grad(U_))) + nonlinearStress_
 274         )
 275     );
 276 }
 277
 278
 279 tmp<fvVectorMatrix> LienCubicKE::divDevReff(volVectorField& U) const
 280 {
 281     return
 282     (
 283         fvc::div(nonlinearStress_)
 284       - fvm::laplacian(nuEff(), U)
 285       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 286     );
 287 }
 288
 289
 290 bool LienCubicKE::read()
 291 {
 292     if (RASModel::read())
 293     {
 294         C1_.readIfPresent(coeffDict());
 295         C2_.readIfPresent(coeffDict());
 296         sigmak_.readIfPresent(coeffDict());
 297         sigmaEps_.readIfPresent(coeffDict());
 298         A1_.readIfPresent(coeffDict());
 299         A2_.readIfPresent(coeffDict());
 300         Ctau1_.readIfPresent(coeffDict());
 301         Ctau2_.readIfPresent(coeffDict());
 302         Ctau3_.readIfPresent(coeffDict());
 303         alphaKsi_.readIfPresent(coeffDict());
 304
 305         return true;
 306     }
 307     else
 308     {
 309         return false;
 310     }
 311 }
 312
 313
 314 void LienCubicKE::correct()
 315 {
 316     RASModel::correct();
 317
 318     if (!turbulence_)
 319     {
 320         return;
 321     }
 322
 323     gradU_ = fvc::grad(U_);
 324
 325     // generation term
 326     volScalarField S2 = symm(gradU_) && gradU_;
 327
 328     volScalarField G
 329     (
 330         "RASModel::G",
 331         Cmu_*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_)
 332     );
 333
 334     // Update espsilon and G at the wall
 335     epsilon_.boundaryField().updateCoeffs();
 336
 337     // Dissipation equation
 338     tmp<fvScalarMatrix> epsEqn
 339     (
 340         fvm::ddt(epsilon_)
 341       + fvm::div(phi_, epsilon_)
 342       - fvm::laplacian(DepsilonEff(), epsilon_)
 343       ==
 344         C1_*G*epsilon_/k_
 345       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 346     );
 347
 348     epsEqn().relax();
 349
 350     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 351
 352     solve(epsEqn);
 353     bound(epsilon_, epsilon0_);
 354
 355
 356     // Turbulent kinetic energy equation
 357
 358     tmp<fvScalarMatrix> kEqn
 359     (
 360         fvm::ddt(k_)
 361       + fvm::div(phi_, k_)
 362       - fvm::laplacian(DkEff(), k_)
 363       ==
 364         G
 365       - fvm::Sp(epsilon_/k_, k_)
 366     );
 367
 368     kEqn().relax();
 369     solve(kEqn);
 370     bound(k_, k0_);
 371
 372
 373     // Re-calculate viscosity
 374
 375     eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())));
 376     ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())));
 377     Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
 378     fEta_ = A2_ + pow(eta_, 3.0);
 379
 380     C5viscosity_ =
 381         - 2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
 382        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()));
 383
 384     nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
 385     nut_.correctBoundaryConditions();
 386
 387     nonlinearStress_ = symm
 388     (
 389         // quadratic terms
 390         pow(k_, 3.0)/sqr(epsilon_)*
 391         (
 392             Ctau1_/fEta_*
 393             (
 394                 (gradU_ & gradU_)
 395               + (gradU_ & gradU_)().T()
 396             )
 397           + Ctau2_/fEta_*(gradU_ & gradU_.T())
 398           + Ctau3_/fEta_*(gradU_.T() & gradU_)
 399         )
 400         // cubic term C4
 401       - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
 402        *pow(Cmu_, 3.0)
 403        *(
 404             ((gradU_ & gradU_) & gradU_.T())
 405           + ((gradU_ & gradU_.T()) & gradU_.T())
 406           - ((gradU_.T() & gradU_) & gradU_)
 407           - ((gradU_.T() & gradU_.T()) & gradU_)
 408         )
 409     );
 410 }
 411
 412
 413 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 414
 415 } // End namespace RASModels
 416 } // End namespace incompressible
 417 } // End namespace Foam
 418
 419 // ************************************************************************* //