src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.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 "kEpsilon.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(kEpsilon, 0);
  43 addToRunTimeSelectionTable(RASModel, kEpsilon, dictionary);
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 kEpsilon::kEpsilon
  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.09
  65         )
  66     ),
  67     C1_
  68     (
  69         dimensioned<scalar>::lookupOrAddToDict
  70         (
  71             "C1",
  72             coeffDict_,
  73             1.44
  74         )
  75     ),
  76     C2_
  77     (
  78         dimensioned<scalar>::lookupOrAddToDict
  79         (
  80             "C2",
  81             coeffDict_,
  82             1.92
  83         )
  84     ),
  85     sigmaEps_
  86     (
  87         dimensioned<scalar>::lookupOrAddToDict
  88         (
  89             "sigmaEps",
  90             coeffDict_,
  91             1.3
  92         )
  93     ),
  94
  95     k_
  96     (
  97         IOobject
  98         (
  99             "k",
 100             runTime_.timeName(),
 101             mesh_,
 102             IOobject::NO_READ,
 103             IOobject::AUTO_WRITE
 104         ),
 105         autoCreateK("k", mesh_)
 106     ),
 107     epsilon_
 108     (
 109         IOobject
 110         (
 111             "epsilon",
 112             runTime_.timeName(),
 113             mesh_,
 114             IOobject::NO_READ,
 115             IOobject::AUTO_WRITE
 116         ),
 117         autoCreateEpsilon("epsilon", mesh_)
 118     ),
 119     nut_
 120     (
 121         IOobject
 122         (
 123             "nut",
 124             runTime_.timeName(),
 125             mesh_,
 126             IOobject::NO_READ,
 127             IOobject::AUTO_WRITE
 128         ),
 129         autoCreateNut("nut", mesh_)
 130     )
 131 {
 132     bound(k_, kMin_);
 133     bound(epsilon_, epsilonMin_);
 134
 135     nut_ = Cmu_*sqr(k_)/epsilon_;
 136     nut_.correctBoundaryConditions();
 137
 138     printCoeffs();
 139 }
 140
 141
 142 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 143
 144 tmp<volSymmTensorField> kEpsilon::R() const
 145 {
 146     return tmp<volSymmTensorField>
 147     (
 148         new volSymmTensorField
 149         (
 150             IOobject
 151             (
 152                 "R",
 153                 runTime_.timeName(),
 154                 mesh_,
 155                 IOobject::NO_READ,
 156                 IOobject::NO_WRITE
 157             ),
 158             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 159             k_.boundaryField().types()
 160         )
 161     );
 162 }
 163
 164
 165 tmp<volSymmTensorField> kEpsilon::devReff() const
 166 {
 167     return tmp<volSymmTensorField>
 168     (
 169         new volSymmTensorField
 170         (
 171             IOobject
 172             (
 173                 "devRhoReff",
 174                 runTime_.timeName(),
 175                 mesh_,
 176                 IOobject::NO_READ,
 177                 IOobject::NO_WRITE
 178             ),
 179            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 180         )
 181     );
 182 }
 183
 184
 185 tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const
 186 {
 187     return
 188     (
 189       - fvm::laplacian(nuEff(), U)
 190       - fvc::div(nuEff()*dev(T(fvc::grad(U))))
 191     );
 192 }
 193
 194
 195 bool kEpsilon::read()
 196 {
 197     if (RASModel::read())
 198     {
 199         Cmu_.readIfPresent(coeffDict());
 200         C1_.readIfPresent(coeffDict());
 201         C2_.readIfPresent(coeffDict());
 202         sigmaEps_.readIfPresent(coeffDict());
 203
 204         return true;
 205     }
 206     else
 207     {
 208         return false;
 209     }
 210 }
 211
 212
 213 void kEpsilon::correct()
 214 {
 215     RASModel::correct();
 216
 217     if (!turbulence_)
 218     {
 219         return;
 220     }
 221
 222     volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_))));
 223
 224     // Update epsilon and G at the wall
 225     epsilon_.boundaryField().updateCoeffs();
 226
 227     // Dissipation equation
 228     tmp<fvScalarMatrix> epsEqn
 229     (
 230         fvm::ddt(epsilon_)
 231       + fvm::div(phi_, epsilon_)
 232       - fvm::Sp(fvc::div(phi_), epsilon_)
 233       - fvm::laplacian(DepsilonEff(), epsilon_)
 234      ==
 235         C1_*G*epsilon_/k_
 236       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 237     );
 238
 239     epsEqn().relax();
 240
 241     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 242
 243     solve(epsEqn);
 244     bound(epsilon_, epsilonMin_);
 245
 246
 247     // Turbulent kinetic energy equation
 248     tmp<fvScalarMatrix> kEqn
 249     (
 250         fvm::ddt(k_)
 251       + fvm::div(phi_, k_)
 252       - fvm::Sp(fvc::div(phi_), k_)
 253       - fvm::laplacian(DkEff(), k_)
 254      ==
 255         G
 256       - fvm::Sp(epsilon_/k_, k_)
 257     );
 258
 259     kEqn().relax();
 260     solve(kEqn);
 261     bound(k_, kMin_);
 262
 263
 264     // Re-calculate viscosity
 265     nut_ = Cmu_*sqr(k_)/epsilon_;
 266     nut_.correctBoundaryConditions();
 267 }
 268
 269
 270 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 271
 272 } // End namespace RASModels
 273 } // End namespace incompressible
 274 } // End namespace Foam
 275
 276 // ************************************************************************* //