src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.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 "realizableKE.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(realizableKE, 0);
  43 addToRunTimeSelectionTable(RASModel, realizableKE, dictionary);
  44
  45 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  46
  47 tmp<volScalarField> realizableKE::rCmu
  48 (
  49     const volTensorField& gradU,
  50     const volScalarField& S2,
  51     const volScalarField& magS
  52 )
  53 {
  54     tmp<volSymmTensorField> tS = dev(symm(gradU));
  55     const volSymmTensorField& S = tS();
  56
  57     volScalarField W =
  58         (2*sqrt(2.0))*((S&S)&&S)
  59        /(
  60             magS*S2
  61           + dimensionedScalar("small", dimensionSet(0, 0, -3, 0, 0), SMALL)
  62         );
  63
  64     tS.clear();
  65
  66     volScalarField phis =
  67         (1.0/3.0)*acos(min(max(sqrt(6.0)*W, -scalar(1)), scalar(1)));
  68     volScalarField As = sqrt(6.0)*cos(phis);
  69     volScalarField Us = sqrt(S2/2.0 + magSqr(skew(gradU)));
  70
  71     return 1.0/(A0_ + As*Us*k_/(epsilon_ + epsilonSmall_));
  72 }
  73
  74
  75 tmp<volScalarField> realizableKE::rCmu
  76 (
  77     const volTensorField& gradU
  78 )
  79 {
  80     volScalarField S2 = 2*magSqr(dev(symm(gradU)));
  81     volScalarField magS = sqrt(S2);
  82     return rCmu(gradU, S2, magS);
  83 }
  84
  85
  86 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  87
  88 realizableKE::realizableKE
  89 (
  90     const volVectorField& U,
  91     const surfaceScalarField& phi,
  92     transportModel& lamTransportModel
  93 )
  94 :
  95     RASModel(typeName, U, phi, lamTransportModel),
  96
  97     Cmu_
  98     (
  99         dimensionedScalar::lookupOrAddToDict
 100         (
 101             "Cmu",
 102             coeffDict_,
 103             0.09
 104         )
 105     ),
 106     A0_
 107     (
 108         dimensionedScalar::lookupOrAddToDict
 109         (
 110             "A0",
 111             coeffDict_,
 112             4.0
 113         )
 114     ),
 115     C2_
 116     (
 117         dimensionedScalar::lookupOrAddToDict
 118         (
 119             "C2",
 120             coeffDict_,
 121             1.9
 122         )
 123     ),
 124     sigmak_
 125     (
 126         dimensionedScalar::lookupOrAddToDict
 127         (
 128             "sigmak",
 129             coeffDict_,
 130             1.0
 131         )
 132     ),
 133     sigmaEps_
 134     (
 135         dimensionedScalar::lookupOrAddToDict
 136         (
 137             "sigmaEps",
 138             coeffDict_,
 139             1.2
 140         )
 141     ),
 142
 143     k_
 144     (
 145         IOobject
 146         (
 147             "k",
 148             runTime_.timeName(),
 149             U_.db(),
 150             IOobject::NO_READ,
 151             IOobject::AUTO_WRITE
 152         ),
 153         autoCreateK("k", mesh_)
 154     ),
 155     epsilon_
 156     (
 157         IOobject
 158         (
 159             "epsilon",
 160             runTime_.timeName(),
 161             U_.db(),
 162             IOobject::NO_READ,
 163             IOobject::AUTO_WRITE
 164         ),
 165         autoCreateEpsilon("epsilon", mesh_)
 166     ),
 167     nut_
 168     (
 169         IOobject
 170         (
 171             "nut",
 172             runTime_.timeName(),
 173             U_.db(),
 174             IOobject::NO_READ,
 175             IOobject::AUTO_WRITE
 176         ),
 177         autoCreateNut("nut", mesh_)
 178     )
 179 {
 180     bound(k_, k0_);
 181     bound(epsilon_, epsilon0_);
 182
 183     nut_ = rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_);
 184     nut_ = min(nut_, nuRatio()*nu());
 185     nut_.correctBoundaryConditions();
 186
 187     printCoeffs();
 188 }
 189
 190
 191 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 192
 193 tmp<volSymmTensorField> realizableKE::R() const
 194 {
 195     return tmp<volSymmTensorField>
 196     (
 197         new volSymmTensorField
 198         (
 199             IOobject
 200             (
 201                 "R",
 202                 runTime_.timeName(),
 203                 U_.db(),
 204                 IOobject::NO_READ,
 205                 IOobject::NO_WRITE
 206             ),
 207             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 208             k_.boundaryField().types()
 209         )
 210     );
 211 }
 212
 213
 214 tmp<volSymmTensorField> realizableKE::devReff() const
 215 {
 216     return tmp<volSymmTensorField>
 217     (
 218         new volSymmTensorField
 219         (
 220             IOobject
 221             (
 222                 "devRhoReff",
 223                 runTime_.timeName(),
 224                 U_.db(),
 225                 IOobject::NO_READ,
 226                 IOobject::NO_WRITE
 227             ),
 228            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 229         )
 230     );
 231 }
 232
 233
 234 tmp<fvVectorMatrix> realizableKE::divDevReff(volVectorField& U) const
 235 {
 236     return
 237     (
 238       - fvm::laplacian(nuEff(), U)
 239       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 240     );
 241 }
 242
 243
 244 bool realizableKE::read()
 245 {
 246     if (RASModel::read())
 247     {
 248         Cmu_.readIfPresent(coeffDict());
 249         A0_.readIfPresent(coeffDict());
 250         C2_.readIfPresent(coeffDict());
 251         sigmak_.readIfPresent(coeffDict());
 252         sigmaEps_.readIfPresent(coeffDict());
 253
 254         return true;
 255     }
 256     else
 257     {
 258         return false;
 259     }
 260 }
 261
 262
 263 void realizableKE::correct()
 264 {
 265     // Bound in case of topological change
 266     // HJ, 22/Aug/2007
 267     if (mesh_.changing())
 268     {
 269         bound(k_, k0_);
 270         bound(epsilon_, epsilon0_);
 271     }
 272
 273     RASModel::correct();
 274
 275     if (!turbulence_)
 276     {
 277         return;
 278     }
 279
 280     volTensorField gradU = fvc::grad(U_);
 281     volScalarField S2 = 2*magSqr(dev(symm(gradU)));
 282     volScalarField magS = sqrt(S2);
 283
 284     volScalarField eta = magS*k_/epsilon_;
 285     volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
 286
 287     volScalarField G("RASModel::G", nut_*S2);
 288
 289     // Update epsilon and G at the wall
 290     epsilon_.boundaryField().updateCoeffs();
 291
 292
 293     // Dissipation equation
 294     tmp<fvScalarMatrix> epsEqn
 295     (
 296         fvm::ddt(epsilon_)
 297       + fvm::div(phi_, epsilon_)
 298       + fvm::SuSp(-fvc::div(phi_), epsilon_)
 299       - fvm::laplacian(DepsilonEff(), epsilon_)
 300      ==
 301         C1*magS*epsilon_
 302       - fvm::Sp
 303         (
 304             C2_*epsilon_/(k_ + sqrt(nu()*epsilon_)),
 305             epsilon_
 306         )
 307     );
 308
 309     epsEqn().relax();
 310
 311     // No longer needed: matrix completes at the point of solution
 312     // HJ, 17/Apr/2012
 313 //     epsEqn().completeAssembly();
 314
 315     solve(epsEqn);
 316     bound(epsilon_, epsilon0_);
 317
 318
 319     // Turbulent kinetic energy equation
 320     tmp<fvScalarMatrix> kEqn
 321     (
 322         fvm::ddt(k_)
 323       + fvm::div(phi_, k_)
 324       + fvm::SuSp(-fvc::div(phi_), k_)
 325       - fvm::laplacian(DkEff(), k_)
 326      ==
 327         G - fvm::Sp(epsilon_/k_, k_)
 328     );
 329
 330     kEqn().relax();
 331     solve(kEqn);
 332     bound(k_, k0_);
 333
 334
 335     // Re-calculate viscosity
 336     nut_ = rCmu(gradU, S2, magS)*sqr(k_)/epsilon_;
 337     nut_ = min(nut_, nuRatio()*nu());
 338     nut_.correctBoundaryConditions();
 339 }
 340
 341
 342 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 343
 344 } // End namespace RASModels
 345 } // End namespace incompressible
 346 } // End namespace Foam
 347
 348 // ************************************************************************* //