src/turbulenceModels/compressible/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 compressible
  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 volScalarField& rho,
  92     const volVectorField& U,
  93     const surfaceScalarField& phi,
  94     const basicThermo& thermophysicalModel
  95 )
  96 :
  97     RASModel(typeName, rho, U, phi, thermophysicalModel),
  98
  99     Cmu_
 100     (
 101         dimensioned<scalar>::lookupOrAddToDict
 102         (
 103             "Cmu",
 104             coeffDict_,
 105             0.09
 106         )
 107     ),
 108     A0_
 109     (
 110         dimensioned<scalar>::lookupOrAddToDict
 111         (
 112             "A0",
 113             coeffDict_,
 114             4.0
 115         )
 116     ),
 117     C2_
 118     (
 119         dimensioned<scalar>::lookupOrAddToDict
 120         (
 121             "C2",
 122             coeffDict_,
 123             1.9
 124         )
 125     ),
 126     sigmak_
 127     (
 128         dimensioned<scalar>::lookupOrAddToDict
 129         (
 130             "sigmak",
 131             coeffDict_,
 132             1.0
 133         )
 134     ),
 135     sigmaEps_
 136     (
 137         dimensioned<scalar>::lookupOrAddToDict
 138         (
 139             "sigmaEps",
 140             coeffDict_,
 141             1.2
 142         )
 143     ),
 144     Prt_
 145     (
 146         dimensioned<scalar>::lookupOrAddToDict
 147         (
 148             "Prt",
 149             coeffDict_,
 150             1.0
 151         )
 152     ),
 153
 154     k_
 155     (
 156         IOobject
 157         (
 158             "k",
 159             runTime_.timeName(),
 160             mesh_,
 161             IOobject::NO_READ,
 162             IOobject::AUTO_WRITE
 163         ),
 164         autoCreateK("k", mesh_)
 165     ),
 166     epsilon_
 167     (
 168         IOobject
 169         (
 170             "epsilon",
 171             runTime_.timeName(),
 172             mesh_,
 173             IOobject::NO_READ,
 174             IOobject::AUTO_WRITE
 175         ),
 176         autoCreateEpsilon("epsilon", mesh_)
 177     ),
 178     mut_
 179     (
 180         IOobject
 181         (
 182             "mut",
 183             runTime_.timeName(),
 184             mesh_,
 185             IOobject::NO_READ,
 186             IOobject::AUTO_WRITE
 187         ),
 188         autoCreateMut("mut", mesh_)
 189     ),
 190     alphat_
 191     (
 192         IOobject
 193         (
 194             "alphat",
 195             runTime_.timeName(),
 196             mesh_,
 197             IOobject::NO_READ,
 198             IOobject::AUTO_WRITE
 199         ),
 200         autoCreateAlphat("alphat", mesh_)
 201     )
 202 {
 203     bound(k_, k0_);
 204     bound(epsilon_, epsilon0_);
 205
 206     mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
 207     mut_.correctBoundaryConditions();
 208
 209     alphat_ = mut_/Prt_;
 210     alphat_.correctBoundaryConditions();
 211
 212     printCoeffs();
 213 }
 214
 215
 216 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 217
 218 tmp<volSymmTensorField> realizableKE::R() const
 219 {
 220     return tmp<volSymmTensorField>
 221     (
 222         new volSymmTensorField
 223         (
 224             IOobject
 225             (
 226                 "R",
 227                 runTime_.timeName(),
 228                 mesh_,
 229                 IOobject::NO_READ,
 230                 IOobject::NO_WRITE
 231             ),
 232             ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
 233             k_.boundaryField().types()
 234         )
 235     );
 236 }
 237
 238
 239 tmp<volSymmTensorField> realizableKE::devRhoReff() const
 240 {
 241     return tmp<volSymmTensorField>
 242     (
 243         new volSymmTensorField
 244         (
 245             IOobject
 246             (
 247                 "devRhoReff",
 248                 runTime_.timeName(),
 249                 mesh_,
 250                 IOobject::NO_READ,
 251                 IOobject::NO_WRITE
 252             ),
 253            -muEff()*dev(twoSymm(fvc::grad(U_)))
 254         )
 255     );
 256 }
 257
 258
 259 tmp<fvVectorMatrix> realizableKE::divDevRhoReff(volVectorField& U) const
 260 {
 261     return
 262     (
 263       - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
 264     );
 265 }
 266
 267
 268 bool realizableKE::read()
 269 {
 270     if (RASModel::read())
 271     {
 272         Cmu_.readIfPresent(coeffDict());
 273         A0_.readIfPresent(coeffDict());
 274         C2_.readIfPresent(coeffDict());
 275         sigmak_.readIfPresent(coeffDict());
 276         sigmaEps_.readIfPresent(coeffDict());
 277         Prt_.readIfPresent(coeffDict());
 278
 279         return true;
 280     }
 281     else
 282     {
 283         return false;
 284     }
 285 }
 286
 287
 288 void realizableKE::correct()
 289 {
 290     if (!turbulence_)
 291     {
 292         // Re-calculate viscosity
 293         mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_;
 294         mut_.correctBoundaryConditions();
 295
 296         // Re-calculate thermal diffusivity
 297         alphat_ = mut_/Prt_;
 298         alphat_.correctBoundaryConditions();
 299
 300         return;
 301     }
 302
 303     RASModel::correct();
 304
 305     volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
 306
 307     if (mesh_.moving())
 308     {
 309         divU += fvc::div(mesh_.phi());
 310     }
 311
 312     volTensorField gradU = fvc::grad(U_);
 313     volScalarField S2 = 2*magSqr(dev(symm(gradU)));
 314     volScalarField magS = sqrt(S2);
 315
 316     volScalarField eta = magS*k_/epsilon_;
 317     volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
 318
 319     volScalarField G("RASModel::G", mut_*(gradU && dev(twoSymm(gradU))));
 320
 321     // Update espsilon and G at the wall
 322     epsilon_.boundaryField().updateCoeffs();
 323
 324     // Dissipation equation
 325     tmp<fvScalarMatrix> epsEqn
 326     (
 327         fvm::ddt(rho_, epsilon_)
 328       + fvm::div(phi_, epsilon_)
 329       - fvm::laplacian(DepsilonEff(), epsilon_)
 330      ==
 331         C1*rho_*magS*epsilon_
 332       - fvm::Sp
 333         (
 334             C2_*rho_*epsilon_/(k_ + sqrt((mu()/rho_)*epsilon_)),
 335             epsilon_
 336         )
 337     );
 338
 339     epsEqn().relax();
 340
 341     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 342
 343     solve(epsEqn);
 344     bound(epsilon_, epsilon0_);
 345
 346
 347     // Turbulent kinetic energy equation
 348
 349     tmp<fvScalarMatrix> kEqn
 350     (
 351         fvm::ddt(rho_, k_)
 352       + fvm::div(phi_, k_)
 353       - fvm::laplacian(DkEff(), k_)
 354      ==
 355         G - fvm::SuSp(2.0/3.0*rho_*divU, k_)
 356       - fvm::Sp(rho_*(epsilon_)/k_, k_)
 357     );
 358
 359     kEqn().relax();
 360     solve(kEqn);
 361     bound(k_, k0_);
 362
 363     // Re-calculate viscosity
 364     mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
 365     mut_.correctBoundaryConditions();
 366
 367     // Re-calculate thermal diffusivity
 368     alphat_ = mut_/Prt_;
 369     alphat_.correctBoundaryConditions();
 370 }
 371
 372
 373 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 374
 375 } // End namespace RASModels
 376 } // End namespace compressible
 377 } // End namespace Foam
 378
 379 // ************************************************************************* //