src/turbulenceModels/RAS/compressible/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 compressible
  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 volScalarField& rho,
  91     const volVectorField& U,
  92     const surfaceScalarField& phi,
  93     basicThermo& thermophysicalModel
  94 )
  95 :
  96     RASModel(typeName, rho, U, phi, thermophysicalModel),
  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     alphak_
 126     (
 127         dimensioned<scalar>::lookupOrAddToDict
 128         (
 129             "alphak",
 130             coeffDict_,
 131             1.0
 132         )
 133     ),
 134     alphaEps_
 135     (
 136         dimensioned<scalar>::lookupOrAddToDict
 137         (
 138             "alphaEps",
 139             coeffDict_,
 140             0.833333
 141         )
 142     ),
 143     alphah_
 144     (
 145         dimensioned<scalar>::lookupOrAddToDict
 146         (
 147             "alphah",
 148             coeffDict_,
 149             1.0
 150         )
 151     ),
 152
 153     k_
 154     (
 155         IOobject
 156         (
 157             "k",
 158             runTime_.timeName(),
 159             mesh_,
 160             IOobject::MUST_READ,
 161             IOobject::AUTO_WRITE
 162         ),
 163         mesh_
 164     ),
 165
 166     epsilon_
 167     (
 168         IOobject
 169         (
 170             "epsilon",
 171             runTime_.timeName(),
 172             mesh_,
 173             IOobject::MUST_READ,
 174             IOobject::AUTO_WRITE
 175         ),
 176         mesh_
 177     ),
 178
 179     mut_
 180     (
 181         IOobject
 182         (
 183             "mut",
 184             runTime_.timeName(),
 185             mesh_,
 186             IOobject::NO_READ,
 187             IOobject::NO_WRITE
 188         ),
 189         rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
 190     )
 191 {
 192     bound(k_, k0_);
 193     bound(epsilon_, epsilon0_);
 194 #   include "wallViscosityI.H"
 195
 196     printCoeffs();
 197 }
 198
 199
 200 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 201
 202 tmp<volSymmTensorField> realizableKE::R() const
 203 {
 204     return tmp<volSymmTensorField>
 205     (
 206         new volSymmTensorField
 207         (
 208             IOobject
 209             (
 210                 "R",
 211                 runTime_.timeName(),
 212                 mesh_,
 213                 IOobject::NO_READ,
 214                 IOobject::NO_WRITE
 215             ),
 216             ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
 217             k_.boundaryField().types()
 218         )
 219     );
 220 }
 221
 222
 223 tmp<volSymmTensorField> realizableKE::devRhoReff() const
 224 {
 225     return tmp<volSymmTensorField>
 226     (
 227         new volSymmTensorField
 228         (
 229             IOobject
 230             (
 231                 "devRhoReff",
 232                 runTime_.timeName(),
 233                 mesh_,
 234                 IOobject::NO_READ,
 235                 IOobject::NO_WRITE
 236             ),
 237            -muEff()*dev(twoSymm(fvc::grad(U_)))
 238         )
 239     );
 240 }
 241
 242
 243 tmp<fvVectorMatrix> realizableKE::divDevRhoReff(volVectorField& U) const
 244 {
 245     return
 246     (
 247       - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
 248     );
 249 }
 250
 251
 252 bool realizableKE::read()
 253 {
 254     if (RASModel::read())
 255     {
 256         Cmu_.readIfPresent(coeffDict_);
 257         A0_.readIfPresent(coeffDict_);
 258         C2_.readIfPresent(coeffDict_);
 259         alphak_.readIfPresent(coeffDict_);
 260         alphaEps_.readIfPresent(coeffDict_);
 261         alphah_.readIfPresent(coeffDict_);
 262
 263         return true;
 264     }
 265     else
 266     {
 267         return false;
 268     }
 269 }
 270
 271
 272 void realizableKE::correct()
 273 {
 274     if (!turbulence_)
 275     {
 276         // Re-calculate viscosity
 277         mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_;
 278         return;
 279     }
 280
 281     RASModel::correct();
 282
 283     volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
 284
 285     if (mesh_.moving())
 286     {
 287         divU += fvc::div(mesh_.phi());
 288     }
 289
 290     volTensorField gradU = fvc::grad(U_);
 291     volScalarField S2 = 2*magSqr(dev(symm(gradU)));
 292     volScalarField magS = sqrt(S2);
 293
 294     volScalarField eta = magS*k_/epsilon_;
 295     volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
 296
 297     volScalarField G = mut_*(gradU && dev(twoSymm(gradU)));
 298
 299 #   include "wallFunctionsI.H"
 300
 301     // Dissipation equation
 302     tmp<fvScalarMatrix> epsEqn
 303     (
 304         fvm::ddt(rho_, epsilon_)
 305       + fvm::div(phi_, epsilon_)
 306       - fvm::laplacian(DepsilonEff(), epsilon_)
 307      ==
 308         C1*rho_*magS*epsilon_
 309       - fvm::Sp
 310         (
 311             C2_*rho_*epsilon_/(k_ + sqrt((mu()/rho_)*epsilon_)),
 312             epsilon_
 313         )
 314     );
 315
 316     epsEqn().relax();
 317
 318 #   include "wallDissipationI.H"
 319
 320     solve(epsEqn);
 321     bound(epsilon_, epsilon0_);
 322
 323
 324     // Turbulent kinetic energy equation
 325
 326     tmp<fvScalarMatrix> kEqn
 327     (
 328         fvm::ddt(rho_, k_)
 329       + fvm::div(phi_, k_)
 330       - fvm::laplacian(DkEff(), k_)
 331      ==
 332         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
 333       - fvm::Sp(rho_*(epsilon_)/k_, k_)
 334     );
 335
 336     kEqn().relax();
 337     solve(kEqn);
 338     bound(k_, k0_);
 339
 340     // Re-calculate viscosity
 341     mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
 342
 343 #   include "wallViscosityI.H"
 344
 345 }
 346
 347
 348 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 349
 350 } // End namespace RASModels
 351 } // End namespace compressible
 352 } // End namespace Foam
 353
 354 // ************************************************************************* //