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