src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.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 "LaunderGibsonRSTM.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(LaunderGibsonRSTM, 0);
  43 addToRunTimeSelectionTable(RASModel, LaunderGibsonRSTM, dictionary);
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 LaunderGibsonRSTM::LaunderGibsonRSTM
  48 (
  49     const volVectorField& U,
  50     const surfaceScalarField& phi,
  51     transportModel& lamTransportModel
  52 )
  53 :
  54     RASModel(typeName, U, phi, lamTransportModel),
  55
  56     Cmu_
  57     (
  58         dimensionedScalar::lookupOrAddToDict
  59         (
  60             "Cmu",
  61             coeffDict_,
  62             0.09
  63         )
  64     ),
  65     kappa_
  66     (
  67         dimensionedScalar::lookupOrAddToDict
  68         (
  69             "kappa",
  70             coeffDict_,
  71             0.41
  72         )
  73     ),
  74     Clg1_
  75     (
  76         dimensionedScalar::lookupOrAddToDict
  77         (
  78             "Clg1",
  79             coeffDict_,
  80             1.8
  81         )
  82     ),
  83     Clg2_
  84     (
  85         dimensionedScalar::lookupOrAddToDict
  86         (
  87             "Clg2",
  88             coeffDict_,
  89             0.6
  90         )
  91     ),
  92     C1_
  93     (
  94         dimensionedScalar::lookupOrAddToDict
  95         (
  96             "C1",
  97             coeffDict_,
  98             1.44
  99         )
 100     ),
 101     C2_
 102     (
 103         dimensionedScalar::lookupOrAddToDict
 104         (
 105             "C2",
 106             coeffDict_,
 107             1.92
 108         )
 109     ),
 110     Cs_
 111     (
 112         dimensionedScalar::lookupOrAddToDict
 113         (
 114             "Cs",
 115             coeffDict_,
 116             0.25
 117         )
 118     ),
 119     Ceps_
 120     (
 121         dimensionedScalar::lookupOrAddToDict
 122         (
 123             "Ceps",
 124             coeffDict_,
 125             0.15
 126         )
 127     ),
 128     sigmaR_
 129     (
 130         dimensionedScalar::lookupOrAddToDict
 131         (
 132             "sigmaR",
 133             coeffDict_,
 134             0.81967
 135         )
 136     ),
 137     sigmaEps_
 138     (
 139         dimensionedScalar::lookupOrAddToDict
 140         (
 141             "sigmaEps",
 142             coeffDict_,
 143             1.3
 144         )
 145     ),
 146     C1Ref_
 147     (
 148         dimensionedScalar::lookupOrAddToDict
 149         (
 150             "C1Ref",
 151             coeffDict_,
 152             0.5
 153         )
 154     ),
 155     C2Ref_
 156     (
 157         dimensionedScalar::lookupOrAddToDict
 158         (
 159             "C2Ref",
 160             coeffDict_,
 161             0.3
 162         )
 163     ),
 164     couplingFactor_
 165     (
 166         dimensionedScalar::lookupOrAddToDict
 167         (
 168             "couplingFactor",
 169             coeffDict_,
 170             0.0
 171         )
 172     ),
 173
 174     yr_(mesh_),
 175
 176     R_
 177     (
 178         IOobject
 179         (
 180             "R",
 181             runTime_.timeName(),
 182             U_.db(),
 183             IOobject::NO_READ,
 184             IOobject::AUTO_WRITE
 185         ),
 186         autoCreateR("R", mesh_)
 187     ),
 188     k_
 189     (
 190         IOobject
 191         (
 192             "k",
 193             runTime_.timeName(),
 194             U_.db(),
 195             IOobject::NO_READ,
 196             IOobject::AUTO_WRITE
 197         ),
 198         autoCreateK("k", mesh_)
 199     ),
 200     epsilon_
 201     (
 202         IOobject
 203         (
 204             "epsilon",
 205             runTime_.timeName(),
 206             U_.db(),
 207             IOobject::NO_READ,
 208             IOobject::AUTO_WRITE
 209         ),
 210         autoCreateEpsilon("epsilon", mesh_)
 211     ),
 212     nut_
 213     (
 214         IOobject
 215         (
 216             "nut",
 217             runTime_.timeName(),
 218             U_.db(),
 219             IOobject::NO_READ,
 220             IOobject::AUTO_WRITE
 221         ),
 222         autoCreateNut("nut", mesh_)
 223     )
 224 {
 225     nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
 226     nut_ = min(nut_, nuRatio()*nu());
 227     nut_.correctBoundaryConditions();
 228
 229     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 230     {
 231         FatalErrorIn
 232         (
 233             "LaunderGibsonRSTM::LaunderGibsonRSTM"
 234             "(const volVectorField& U, const surfaceScalarField& phi,"
 235             "transportModel& lamTransportModel)"
 236         )   << "couplingFactor = " << couplingFactor_
 237             << " is not in range 0 - 1" << nl
 238             << exit(FatalError);
 239     }
 240
 241     printCoeffs();
 242 }
 243
 244
 245 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 246
 247 tmp<volSymmTensorField> LaunderGibsonRSTM::devReff() const
 248 {
 249     return tmp<volSymmTensorField>
 250     (
 251         new volSymmTensorField
 252         (
 253             IOobject
 254             (
 255                 "devRhoReff",
 256                 runTime_.timeName(),
 257                 U_.db(),
 258                 IOobject::NO_READ,
 259                 IOobject::NO_WRITE
 260             ),
 261             R_ - nu()*dev(twoSymm(fvc::grad(U_)))
 262         )
 263     );
 264 }
 265
 266
 267 tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
 268 {
 269     if (couplingFactor_.value() > 0.0)
 270     {
 271         return
 272         (
 273             fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
 274           + fvc::laplacian((1.0-couplingFactor_)*nut_, U, "laplacian(nuEff,U)")
 275           - fvm::laplacian(nuEff(), U)
 276         );
 277     }
 278     else
 279     {
 280         return
 281         (
 282             fvc::div(R_)
 283           + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
 284           - fvm::laplacian(nuEff(), U)
 285         );
 286     }
 287 }
 288
 289
 290 bool LaunderGibsonRSTM::read()
 291 {
 292     if (RASModel::read())
 293     {
 294         Cmu_.readIfPresent(coeffDict());
 295         Clg1_.readIfPresent(coeffDict());
 296         Clg2_.readIfPresent(coeffDict());
 297         C1_.readIfPresent(coeffDict());
 298         C2_.readIfPresent(coeffDict());
 299         Cs_.readIfPresent(coeffDict());
 300         Ceps_.readIfPresent(coeffDict());
 301         sigmaR_.readIfPresent(coeffDict());
 302         sigmaEps_.readIfPresent(coeffDict());
 303         C1Ref_.readIfPresent(coeffDict());
 304         C2Ref_.readIfPresent(coeffDict());
 305
 306         couplingFactor_.readIfPresent(coeffDict());
 307
 308         if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 309         {
 310             FatalErrorIn("LaunderGibsonRSTM::read()")
 311                 << "couplingFactor = " << couplingFactor_
 312                 << " is not in range 0 - 1"
 313                 << exit(FatalError);
 314         }
 315
 316         return true;
 317     }
 318     else
 319     {
 320         return false;
 321     }
 322 }
 323
 324
 325 void LaunderGibsonRSTM::correct()
 326 {
 327     // Bound in case of topological change
 328     // HJ, 22/Aug/2007
 329     if (mesh_.changing())
 330     {
 331         R_.correctBoundaryConditions();
 332         bound(epsilon_, epsilon0_);
 333     }
 334
 335     RASModel::correct();
 336
 337     if (!turbulence_)
 338     {
 339         return;
 340     }
 341
 342     if (mesh_.changing())
 343     {
 344         yr_.correct();
 345     }
 346
 347     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
 348     volSymmTensorField C = -fvc::div(phi_, R_);
 349     volScalarField G("RASModel::G", 0.5*mag(tr(P)));
 350
 351     // Update epsilon and G at the wall
 352     epsilon_.boundaryField().updateCoeffs();
 353
 354     // Dissipation equation
 355     tmp<fvScalarMatrix> epsEqn
 356     (
 357         fvm::ddt(epsilon_)
 358       + fvm::div(phi_, epsilon_)
 359       + fvm::SuSp(-fvc::div(phi_), epsilon_)
 360     //- fvm::laplacian(Ceps*(k_/epsilon_)*R_, epsilon_)
 361       - fvm::laplacian(DepsilonEff(), epsilon_)
 362      ==
 363         C1_*G*epsilon_/k_
 364       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 365     );
 366
 367     epsEqn().relax();
 368
 369     // No longer needed: matrix completes at the point of solution
 370     // HJ, 17/Apr/2012
 371 //     epsEqn().completeAssembly();
 372
 373     solve(epsEqn);
 374     bound(epsilon_, epsilon0_);
 375
 376
 377     // Reynolds stress equation
 378
 379     const fvPatchList& patches = mesh_.boundary();
 380
 381     forAll(patches, patchi)
 382     {
 383         const fvPatch& curPatch = patches[patchi];
 384
 385         if (curPatch.isWall())
 386         {
 387             forAll(curPatch, facei)
 388             {
 389                 label faceCelli = curPatch.faceCells()[facei];
 390                 P[faceCelli] *=
 391                     min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 1.0);
 392             }
 393         }
 394     }
 395
 396     volSymmTensorField reflect = C1Ref_*epsilon_/k_*R_ - C2Ref_*Clg2_*dev(P);
 397
 398     tmp<fvSymmTensorMatrix> REqn
 399     (
 400         fvm::ddt(R_)
 401       + fvm::div(phi_, R_)
 402       + fvm::SuSp(-fvc::div(phi_), R_)
 403     //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
 404       - fvm::laplacian(DREff(), R_)
 405       + fvm::Sp(Clg1_*epsilon_/k_, R_)
 406       ==
 407         P
 408       + (2.0/3.0*(Clg1_ - 1)*I)*epsilon_
 409         // Change for consistency with Fluent implementation.
 410         // Emil Baric, NUMAP-FOAM 2011
 411         // HJ, 13/Dec/2011
 412       - Clg2_*(dev(P) - dev(C))
 413
 414         // wall reflection terms
 415       + symm
 416         (
 417             I*((yr_.n() & reflect) & yr_.n())
 418           - 1.5*(yr_.n()*(reflect & yr_.n())
 419           + (yr_.n() & reflect)*yr_.n())
 420         )*pow(Cmu_, 0.75)*pow(k_, 1.5)/(kappa_*yr_*epsilon_)
 421     );
 422
 423     REqn().relax();
 424     solve(REqn);
 425
 426     R_.max
 427     (
 428         dimensionedSymmTensor
 429         (
 430             "zero",
 431             R_.dimensions(),
 432             symmTensor
 433             (
 434                 k0_.value(), -GREAT, -GREAT,
 435                              k0_.value(), -GREAT,
 436                                           k0_.value()
 437             )
 438         )
 439     );
 440
 441     k_ == 0.5*tr(R_);
 442     bound(k_, k0_);
 443
 444
 445     // Re-calculate turbulent viscosity
 446     nut_ = Cmu_*sqr(k_)/epsilon_;
 447     nut_ = min(nut_, nuRatio()*nu());
 448     nut_.correctBoundaryConditions();
 449
 450
 451     // Correct wall shear stresses.  Moved to wall functions with anisotropy
 452     // HJ, 14/Dec/2011
 453 }
 454
 455
 456 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 457
 458 } // End namespace RASModels
 459 } // End namespace incompressible
 460 } // End namespace Foam
 461
 462 // ************************************************************************* //