src/turbulenceModels/RAS/incompressible/kEpsilon/kEpsilon.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "kEpsilon.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "wallFvPatch.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 // Construct from components
  48 kEpsilon::kEpsilon
  49 (
  50     const volVectorField& U,
  51     const surfaceScalarField& phi,
  52     transportModel& lamTransportModel
  53 )
  54 :
  55     RASModel(typeName, U, phi, lamTransportModel),
  56
  57     Cmu_
  58     (
  59         dimensioned<scalar>::lookupOrAddToDict
  60         (
  61             "Cmu",
  62             coeffDict_,
  63             0.09
  64         )
  65     ),
  66     C1_
  67     (
  68         dimensioned<scalar>::lookupOrAddToDict
  69         (
  70             "C1",
  71             coeffDict_,
  72             1.44
  73         )
  74     ),
  75     C2_
  76     (
  77         dimensioned<scalar>::lookupOrAddToDict
  78         (
  79             "C2",
  80             coeffDict_,
  81             1.92
  82         )
  83     ),
  84     alphaEps_
  85     (
  86         dimensioned<scalar>::lookupOrAddToDict
  87         (
  88             "alphaEps",
  89             coeffDict_,
  90             0.76923
  91         )
  92     ),
  93
  94     k_
  95     (
  96         IOobject
  97         (
  98             "k",
  99             runTime_.timeName(),
 100             mesh_,
 101             IOobject::MUST_READ,
 102             IOobject::AUTO_WRITE
 103         ),
 104         mesh_
 105     ),
 106
 107     epsilon_
 108     (
 109         IOobject
 110         (
 111             "epsilon",
 112             runTime_.timeName(),
 113             mesh_,
 114             IOobject::MUST_READ,
 115             IOobject::AUTO_WRITE
 116         ),
 117         mesh_
 118     ),
 119
 120     nut_(Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_))
 121 {
 122 #   include "wallViscosityI.H"
 123
 124     printCoeffs();
 125 }
 126
 127
 128 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 129
 130 tmp<volSymmTensorField> kEpsilon::R() const
 131 {
 132     return tmp<volSymmTensorField>
 133     (
 134         new volSymmTensorField
 135         (
 136             IOobject
 137             (
 138                 "R",
 139                 runTime_.timeName(),
 140                 mesh_,
 141                 IOobject::NO_READ,
 142                 IOobject::NO_WRITE
 143             ),
 144             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 145             k_.boundaryField().types()
 146         )
 147     );
 148 }
 149
 150
 151 tmp<volSymmTensorField> kEpsilon::devReff() const
 152 {
 153     return tmp<volSymmTensorField>
 154     (
 155         new volSymmTensorField
 156         (
 157             IOobject
 158             (
 159                 "devRhoReff",
 160                 runTime_.timeName(),
 161                 mesh_,
 162                 IOobject::NO_READ,
 163                 IOobject::NO_WRITE
 164             ),
 165            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 166         )
 167     );
 168 }
 169
 170
 171 tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const
 172 {
 173     return
 174     (
 175       - fvm::laplacian(nuEff(), U)
 176       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 177     );
 178 }
 179
 180
 181 bool kEpsilon::read()
 182 {
 183     if (RASModel::read())
 184     {
 185         Cmu_.readIfPresent(coeffDict_);
 186         C1_.readIfPresent(coeffDict_);
 187         C2_.readIfPresent(coeffDict_);
 188         alphaEps_.readIfPresent(coeffDict_);
 189
 190         return true;
 191     }
 192     else
 193     {
 194         return false;
 195     }
 196 }
 197
 198
 199 void kEpsilon::correct()
 200 {
 201     transportModel_.correct();
 202
 203     if (!turbulence_)
 204     {
 205         return;
 206     }
 207
 208     RASModel::correct();
 209
 210     volScalarField G = nut_*2*magSqr(symm(fvc::grad(U_)));
 211
 212 #   include "wallFunctionsI.H"
 213
 214     // Dissipation equation
 215     tmp<fvScalarMatrix> epsEqn
 216     (
 217         fvm::ddt(epsilon_)
 218       + fvm::div(phi_, epsilon_)
 219       - fvm::Sp(fvc::div(phi_), epsilon_)
 220       - fvm::laplacian(DepsilonEff(), epsilon_)
 221      ==
 222         C1_*G*epsilon_/k_
 223       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 224     );
 225
 226     epsEqn().relax();
 227
 228 #   include "wallDissipationI.H"
 229
 230     solve(epsEqn);
 231     bound(epsilon_, epsilon0_);
 232
 233
 234     // Turbulent kinetic energy equation
 235     tmp<fvScalarMatrix> kEqn
 236     (
 237         fvm::ddt(k_)
 238       + fvm::div(phi_, k_)
 239       - fvm::Sp(fvc::div(phi_), k_)
 240       - fvm::laplacian(DkEff(), k_)
 241      ==
 242         G
 243       - fvm::Sp(epsilon_/k_, k_)
 244     );
 245
 246     kEqn().relax();
 247     solve(kEqn);
 248     bound(k_, k0_);
 249
 250
 251     // Re-calculate viscosity
 252     nut_ = Cmu_*sqr(k_)/epsilon_;
 253
 254 #   include "wallViscosityI.H"
 255
 256 }
 257
 258
 259 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 260
 261 } // End namespace RASModels
 262 } // End namespace incompressible
 263 } // End namespace Foam
 264
 265 // ************************************************************************* //