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