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