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