src/turbulenceModels/RAS/incompressible/realizableKE/realizableKE.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 "realizableKE.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(realizableKE, 0);
  43 addToRunTimeSelectionTable(RASModel, realizableKE, dictionary);
  44
  45 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  46
  47 tmp<volScalarField> realizableKE::rCmu
  48 (
  49     const volTensorField& gradU,
  50     const volScalarField& S2,
  51     const volScalarField& magS
  52 )
  53 {
  54     tmp<volSymmTensorField> tS = dev(symm(gradU));
  55     const volSymmTensorField& S = tS();
  56
  57     volScalarField W =
  58         (2*sqrt(2.0))*((S&S)&&S)
  59        /(
  60             magS*S2
  61           + dimensionedScalar("small", dimensionSet(0, 0, -3, 0, 0), SMALL)
  62         );
  63
  64     tS.clear();
  65
  66     volScalarField phis =
  67         (1.0/3.0)*acos(min(max(sqrt(6.0)*W, -scalar(1)), scalar(1)));
  68     volScalarField As = sqrt(6.0)*cos(phis);
  69     volScalarField Us = sqrt(S2/2.0 + magSqr(skew(gradU)));
  70
  71     return 1.0/(A0_ + As*Us*k_/(epsilon_ + epsilonSmall_));
  72 }
  73
  74
  75 tmp<volScalarField> realizableKE::rCmu
  76 (
  77     const volTensorField& gradU
  78 )
  79 {
  80     volScalarField S2 = 2*magSqr(dev(symm(gradU)));
  81     volScalarField magS = sqrt(S2);
  82     return rCmu(gradU, S2, magS);
  83 }
  84
  85
  86 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  87
  88 realizableKE::realizableKE
  89 (
  90     const volVectorField& U,
  91     const surfaceScalarField& phi,
  92     transportModel& lamTransportModel
  93 )
  94 :
  95     RASModel(typeName, U, phi, lamTransportModel),
  96
  97     Cmu_
  98     (
  99         dimensioned<scalar>::lookupOrAddToDict
 100         (
 101             "Cmu",
 102             coeffDict_,
 103             0.09
 104         )
 105     ),
 106     A0_
 107     (
 108         dimensioned<scalar>::lookupOrAddToDict
 109         (
 110             "A0",
 111             coeffDict_,
 112             4.0
 113         )
 114     ),
 115     C2_
 116     (
 117         dimensioned<scalar>::lookupOrAddToDict
 118         (
 119             "C2",
 120             coeffDict_,
 121             1.9
 122         )
 123     ),
 124     alphak_
 125     (
 126         dimensioned<scalar>::lookupOrAddToDict
 127         (
 128             "alphak",
 129             coeffDict_,
 130             1.0
 131         )
 132     ),
 133     alphaEps_
 134     (
 135         dimensioned<scalar>::lookupOrAddToDict
 136         (
 137             "alphaEps",
 138             coeffDict_,
 139             0.833333
 140         )
 141     ),
 142
 143     k_
 144     (
 145         IOobject
 146         (
 147             "k",
 148             runTime_.timeName(),
 149             mesh_,
 150             IOobject::MUST_READ,
 151             IOobject::AUTO_WRITE
 152         ),
 153         mesh_
 154     ),
 155
 156     epsilon_
 157     (
 158         IOobject
 159         (
 160             "epsilon",
 161             runTime_.timeName(),
 162             mesh_,
 163             IOobject::MUST_READ,
 164             IOobject::AUTO_WRITE
 165         ),
 166         mesh_
 167     ),
 168
 169     nut_(rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_))
 170 {
 171     bound(k_, k0_);
 172     bound(epsilon_, epsilon0_);
 173 #   include "wallViscosityI.H"
 174
 175     printCoeffs();
 176 }
 177
 178
 179 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 180
 181 tmp<volSymmTensorField> realizableKE::R() const
 182 {
 183     return tmp<volSymmTensorField>
 184     (
 185         new volSymmTensorField
 186         (
 187             IOobject
 188             (
 189                 "R",
 190                 runTime_.timeName(),
 191                 mesh_,
 192                 IOobject::NO_READ,
 193                 IOobject::NO_WRITE
 194             ),
 195             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 196             k_.boundaryField().types()
 197         )
 198     );
 199 }
 200
 201
 202 tmp<volSymmTensorField> realizableKE::devReff() const
 203 {
 204     return tmp<volSymmTensorField>
 205     (
 206         new volSymmTensorField
 207         (
 208             IOobject
 209             (
 210                 "devRhoReff",
 211                 runTime_.timeName(),
 212                 mesh_,
 213                 IOobject::NO_READ,
 214                 IOobject::NO_WRITE
 215             ),
 216            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 217         )
 218     );
 219 }
 220
 221
 222 tmp<fvVectorMatrix> realizableKE::divDevReff(volVectorField& U) const
 223 {
 224     return
 225     (
 226       - fvm::laplacian(nuEff(), U)
 227       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 228     );
 229 }
 230
 231
 232 bool realizableKE::read()
 233 {
 234     if (RASModel::read())
 235     {
 236         Cmu_.readIfPresent(coeffDict_);
 237         A0_.readIfPresent(coeffDict_);
 238         C2_.readIfPresent(coeffDict_);
 239         alphak_.readIfPresent(coeffDict_);
 240         alphaEps_.readIfPresent(coeffDict_);
 241
 242         return true;
 243     }
 244     else
 245     {
 246         return false;
 247     }
 248 }
 249
 250
 251 void realizableKE::correct()
 252 {
 253     transportModel_.correct();
 254
 255     if (!turbulence_)
 256     {
 257         return;
 258     }
 259
 260     RASModel::correct();
 261
 262     volTensorField gradU = fvc::grad(U_);
 263     volScalarField S2 = 2*magSqr(dev(symm(gradU)));
 264     volScalarField magS = sqrt(S2);
 265
 266     volScalarField eta = magS*k_/epsilon_;
 267     volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
 268
 269     volScalarField G = nut_*S2;
 270
 271 #   include "wallFunctionsI.H"
 272
 273
 274     // Dissipation equation
 275     tmp<fvScalarMatrix> epsEqn
 276     (
 277         fvm::ddt(epsilon_)
 278       + fvm::div(phi_, epsilon_)
 279       - fvm::Sp(fvc::div(phi_), epsilon_)
 280       - fvm::laplacian(DepsilonEff(), epsilon_)
 281      ==
 282         C1*magS*epsilon_
 283       - fvm::Sp
 284         (
 285             C2_*epsilon_/(k_ + sqrt(nu()*epsilon_)),
 286             epsilon_
 287         )
 288     );
 289
 290     epsEqn().relax();
 291
 292 #   include "wallDissipationI.H"
 293
 294     solve(epsEqn);
 295     bound(epsilon_, epsilon0_);
 296
 297
 298     // Turbulent kinetic energy equation
 299     tmp<fvScalarMatrix> kEqn
 300     (
 301         fvm::ddt(k_)
 302       + fvm::div(phi_, k_)
 303       - fvm::Sp(fvc::div(phi_), k_)
 304       - fvm::laplacian(DkEff(), k_)
 305      ==
 306         G - fvm::Sp(epsilon_/k_, k_)
 307     );
 308
 309     kEqn().relax();
 310     solve(kEqn);
 311     bound(k_, k0_);
 312
 313
 314     // Re-calculate viscosity
 315     nut_ = rCmu(gradU, S2, magS)*sqr(k_)/epsilon_;
 316
 317 #   include "wallViscosityI.H"
 318
 319 }
 320
 321
 322 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 323
 324 } // End namespace RASModels
 325 } // End namespace incompressible
 326 } // End namespace Foam
 327
 328 // ************************************************************************* //