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