src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, 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& transport,
  52     const word& turbulenceModelName,
  53     const word& modelName
  54 )
  55 :
  56     RASModel(modelName, U, phi, transport, turbulenceModelName),
  57
  58     Cmu_
  59     (
  60         dimensioned<scalar>::lookupOrAddToDict
  61         (
  62             "Cmu",
  63             coeffDict_,
  64             0.0845
  65         )
  66     ),
  67     C1_
  68     (
  69         dimensioned<scalar>::lookupOrAddToDict
  70         (
  71             "C1",
  72             coeffDict_,
  73             1.42
  74         )
  75     ),
  76     C2_
  77     (
  78         dimensioned<scalar>::lookupOrAddToDict
  79         (
  80             "C2",
  81             coeffDict_,
  82             1.68
  83         )
  84     ),
  85     sigmak_
  86     (
  87         dimensioned<scalar>::lookupOrAddToDict
  88         (
  89             "sigmak",
  90             coeffDict_,
  91             0.71942
  92         )
  93     ),
  94     sigmaEps_
  95     (
  96         dimensioned<scalar>::lookupOrAddToDict
  97         (
  98             "sigmaEps",
  99             coeffDict_,
 100             0.71942
 101         )
 102     ),
 103     eta0_
 104     (
 105         dimensioned<scalar>::lookupOrAddToDict
 106         (
 107             "eta0",
 108             coeffDict_,
 109             4.38
 110         )
 111     ),
 112     beta_
 113     (
 114         dimensioned<scalar>::lookupOrAddToDict
 115         (
 116             "beta",
 117             coeffDict_,
 118             0.012
 119         )
 120     ),
 121
 122     k_
 123     (
 124         IOobject
 125         (
 126             "k",
 127             runTime_.timeName(),
 128             mesh_,
 129             IOobject::NO_READ,
 130             IOobject::AUTO_WRITE
 131         ),
 132         autoCreateK("k", mesh_)
 133     ),
 134     epsilon_
 135     (
 136         IOobject
 137         (
 138             "epsilon",
 139             runTime_.timeName(),
 140             mesh_,
 141             IOobject::NO_READ,
 142             IOobject::AUTO_WRITE
 143         ),
 144         autoCreateEpsilon("epsilon", mesh_)
 145     ),
 146     nut_
 147     (
 148         IOobject
 149         (
 150             "nut",
 151             runTime_.timeName(),
 152             mesh_,
 153             IOobject::NO_READ,
 154             IOobject::AUTO_WRITE
 155         ),
 156         autoCreateNut("nut", mesh_)
 157     )
 158 {
 159     bound(k_, kMin_);
 160     bound(epsilon_, epsilonMin_);
 161
 162     nut_ = Cmu_*sqr(k_)/epsilon_;
 163     nut_.correctBoundaryConditions();
 164
 165     printCoeffs();
 166 }
 167
 168
 169 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 170
 171
 172 tmp<volSymmTensorField> RNGkEpsilon::R() const
 173 {
 174     return tmp<volSymmTensorField>
 175     (
 176         new volSymmTensorField
 177         (
 178             IOobject
 179             (
 180                 "R",
 181                 runTime_.timeName(),
 182                 mesh_,
 183                 IOobject::NO_READ,
 184                 IOobject::NO_WRITE
 185             ),
 186             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 187             k_.boundaryField().types()
 188         )
 189     );
 190 }
 191
 192
 193 tmp<volSymmTensorField> RNGkEpsilon::devReff() const
 194 {
 195     return tmp<volSymmTensorField>
 196     (
 197         new volSymmTensorField
 198         (
 199             IOobject
 200             (
 201                 "devRhoReff",
 202                 runTime_.timeName(),
 203                 mesh_,
 204                 IOobject::NO_READ,
 205                 IOobject::NO_WRITE
 206             ),
 207            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 208         )
 209     );
 210 }
 211
 212
 213 tmp<fvVectorMatrix> RNGkEpsilon::divDevReff(volVectorField& U) const
 214 {
 215     return
 216     (
 217       - fvm::laplacian(nuEff(), U)
 218       - fvc::div(nuEff()*dev(T(fvc::grad(U))))
 219     );
 220 }
 221
 222
 223 bool RNGkEpsilon::read()
 224 {
 225     if (RASModel::read())
 226     {
 227         Cmu_.readIfPresent(coeffDict());
 228         C1_.readIfPresent(coeffDict());
 229         C2_.readIfPresent(coeffDict());
 230         sigmak_.readIfPresent(coeffDict());
 231         sigmaEps_.readIfPresent(coeffDict());
 232         eta0_.readIfPresent(coeffDict());
 233         beta_.readIfPresent(coeffDict());
 234
 235         return true;
 236     }
 237     else
 238     {
 239         return false;
 240     }
 241 }
 242
 243
 244 void RNGkEpsilon::correct()
 245 {
 246     RASModel::correct();
 247
 248     if (!turbulence_)
 249     {
 250         return;
 251     }
 252
 253     const volScalarField S2(2*magSqr(symm(fvc::grad(U_))));
 254     volScalarField G("RASModel::G", nut_*S2);
 255
 256     const volScalarField eta(sqrt(S2)*k_/epsilon_);
 257     volScalarField R
 258     (
 259         ((eta*(scalar(1) - eta/eta0_))/(scalar(1) + beta_*eta*sqr(eta)))
 260     );
 261
 262     // Update epsilon and G at the wall
 263     epsilon_.boundaryField().updateCoeffs();
 264
 265     // Dissipation equation
 266     tmp<fvScalarMatrix> epsEqn
 267     (
 268         fvm::ddt(epsilon_)
 269       + fvm::div(phi_, epsilon_)
 270       - fvm::laplacian(DepsilonEff(), epsilon_)
 271      ==
 272         (C1_ - R)*G*epsilon_/k_
 273     //- fvm::SuSp(R*G/k_, epsilon_)
 274       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 275     );
 276
 277     epsEqn().relax();
 278
 279     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 280
 281     solve(epsEqn);
 282     bound(epsilon_, epsilonMin_);
 283
 284
 285     // Turbulent kinetic energy equation
 286
 287     tmp<fvScalarMatrix> kEqn
 288     (
 289         fvm::ddt(k_)
 290       + fvm::div(phi_, k_)
 291       - fvm::laplacian(DkEff(), k_)
 292      ==
 293         G - fvm::Sp(epsilon_/k_, k_)
 294     );
 295
 296     kEqn().relax();
 297     solve(kEqn);
 298     bound(k_, kMin_);
 299
 300
 301     // Re-calculate viscosity
 302     nut_ = Cmu_*sqr(k_)/epsilon_;
 303     nut_.correctBoundaryConditions();
 304 }
 305
 306
 307 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 308
 309 } // End namespace RASModels
 310 } // End namespace incompressible
 311 } // End namespace Foam
 312
 313 // ************************************************************************* //