src/turbulenceModels/compressible/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 compressible
  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 volScalarField& rho,
  51     const volVectorField& U,
  52     const surfaceScalarField& phi,
  53     const basicThermo& thermophysicalModel
  54 )
  55 :
  56     RASModel(typeName, rho, U, phi, thermophysicalModel),
  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     C3_
  86     (
  87         dimensioned<scalar>::lookupOrAddToDict
  88         (
  89             "C3",
  90             coeffDict_,
  91             -0.33
  92         )
  93     ),
  94     sigmak_
  95     (
  96         dimensioned<scalar>::lookupOrAddToDict
  97         (
  98             "sigmak",
  99             coeffDict_,
 100             1.0
 101         )
 102     ),
 103     sigmaEps_
 104     (
 105         dimensioned<scalar>::lookupOrAddToDict
 106         (
 107             "sigmaEps",
 108             coeffDict_,
 109             1.3
 110         )
 111     ),
 112     Prt_
 113     (
 114         dimensioned<scalar>::lookupOrAddToDict
 115         (
 116             "Prt",
 117             coeffDict_,
 118             1.0
 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     mut_
 147     (
 148         IOobject
 149         (
 150             "mut",
 151             runTime_.timeName(),
 152             mesh_,
 153             IOobject::NO_READ,
 154             IOobject::AUTO_WRITE
 155         ),
 156         autoCreateMut("mut", mesh_)
 157     ),
 158     alphat_
 159     (
 160         IOobject
 161         (
 162             "alphat",
 163             runTime_.timeName(),
 164             mesh_,
 165             IOobject::NO_READ,
 166             IOobject::AUTO_WRITE
 167         ),
 168         autoCreateAlphat("alphat", mesh_)
 169     )
 170 {
 171     mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
 172     mut_.correctBoundaryConditions();
 173
 174     alphat_ = mut_/Prt_;
 175     alphat_.correctBoundaryConditions();
 176
 177     printCoeffs();
 178 }
 179
 180
 181 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 182
 183 tmp<volSymmTensorField> kEpsilon::R() const
 184 {
 185     return tmp<volSymmTensorField>
 186     (
 187         new volSymmTensorField
 188         (
 189             IOobject
 190             (
 191                 "R",
 192                 runTime_.timeName(),
 193                 mesh_,
 194                 IOobject::NO_READ,
 195                 IOobject::NO_WRITE
 196             ),
 197             ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
 198             k_.boundaryField().types()
 199         )
 200     );
 201 }
 202
 203
 204 tmp<volSymmTensorField> kEpsilon::devRhoReff() const
 205 {
 206     return tmp<volSymmTensorField>
 207     (
 208         new volSymmTensorField
 209         (
 210             IOobject
 211             (
 212                 "devRhoReff",
 213                 runTime_.timeName(),
 214                 mesh_,
 215                 IOobject::NO_READ,
 216                 IOobject::NO_WRITE
 217             ),
 218            -muEff()*dev(twoSymm(fvc::grad(U_)))
 219         )
 220     );
 221 }
 222
 223
 224 tmp<fvVectorMatrix> kEpsilon::divDevRhoReff(volVectorField& U) const
 225 {
 226     return
 227     (
 228       - fvm::laplacian(muEff(), U)
 229       - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
 230     );
 231 }
 232
 233
 234 bool kEpsilon::read()
 235 {
 236     if (RASModel::read())
 237     {
 238         Cmu_.readIfPresent(coeffDict());
 239         C1_.readIfPresent(coeffDict());
 240         C2_.readIfPresent(coeffDict());
 241         C3_.readIfPresent(coeffDict());
 242         sigmak_.readIfPresent(coeffDict());
 243         sigmaEps_.readIfPresent(coeffDict());
 244         Prt_.readIfPresent(coeffDict());
 245
 246         return true;
 247     }
 248     else
 249     {
 250         return false;
 251     }
 252 }
 253
 254
 255 void kEpsilon::correct()
 256 {
 257     if (!turbulence_)
 258     {
 259         // Re-calculate viscosity
 260         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
 261         mut_.correctBoundaryConditions();
 262
 263         // Re-calculate thermal diffusivity
 264         alphat_ = mut_/Prt_;
 265         alphat_.correctBoundaryConditions();
 266
 267         return;
 268     }
 269
 270     RASModel::correct();
 271
 272     volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
 273
 274     if (mesh_.moving())
 275     {
 276         divU += fvc::div(mesh_.phi());
 277     }
 278
 279     tmp<volTensorField> tgradU = fvc::grad(U_);
 280     volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
 281     tgradU.clear();
 282
 283     // Update espsilon and G at the wall
 284     epsilon_.boundaryField().updateCoeffs();
 285
 286     // Dissipation equation
 287     tmp<fvScalarMatrix> epsEqn
 288     (
 289         fvm::ddt(rho_, epsilon_)
 290       + fvm::div(phi_, epsilon_)
 291       - fvm::laplacian(DepsilonEff(), epsilon_)
 292      ==
 293         C1_*G*epsilon_/k_
 294       - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)
 295       - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
 296     );
 297
 298     epsEqn().relax();
 299
 300     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 301
 302     solve(epsEqn);
 303     bound(epsilon_, epsilon0_);
 304
 305
 306     // Turbulent kinetic energy equation
 307
 308     tmp<fvScalarMatrix> kEqn
 309     (
 310         fvm::ddt(rho_, k_)
 311       + fvm::div(phi_, k_)
 312       - fvm::laplacian(DkEff(), k_)
 313      ==
 314         G
 315       - fvm::SuSp((2.0/3.0)*rho_*divU, k_)
 316       - fvm::Sp(rho_*epsilon_/k_, k_)
 317     );
 318
 319     kEqn().relax();
 320     solve(kEqn);
 321     bound(k_, k0_);
 322
 323
 324     // Re-calculate viscosity
 325     mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
 326     mut_.correctBoundaryConditions();
 327
 328     // Re-calculate thermal diffusivity
 329     alphat_ = mut_/Prt_;
 330     alphat_.correctBoundaryConditions();
 331 }
 332
 333
 334 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 335
 336 } // End namespace RASModels
 337 } // End namespace compressible
 338 } // End namespace Foam
 339
 340 // ************************************************************************* //