src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  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 "SpalartAllmaras.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(SpalartAllmaras, 0);
  43 addToRunTimeSelectionTable(RASModel, SpalartAllmaras, dictionary);
  44
  45 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  46
  47 tmp<volScalarField> SpalartAllmaras::chi() const
  48 {
  49     return rho_*nuTilda_/mu();
  50 }
  51
  52
  53 tmp<volScalarField> SpalartAllmaras::fv1(const volScalarField& chi) const
  54 {
  55     volScalarField chi3(pow3(chi));
  56     return chi3/(chi3 + pow3(Cv1_));
  57 }
  58
  59
  60 tmp<volScalarField> SpalartAllmaras::fv2
  61 (
  62     const volScalarField& chi,
  63     const volScalarField& fv1
  64 ) const
  65 {
  66     return 1.0/pow3(scalar(1) + chi/Cv2_);
  67 }
  68
  69
  70 tmp<volScalarField> SpalartAllmaras::fv3
  71 (
  72     const volScalarField& chi,
  73     const volScalarField& fv1
  74 ) const
  75 {
  76     volScalarField chiByCv2((1/Cv2_)*chi);
  77
  78     return
  79         (scalar(1) + chi*fv1)
  80        *(1/Cv2_)
  81        *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
  82        /pow3(scalar(1) + chiByCv2);
  83 }
  84
  85
  86 tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& Stilda) const
  87 {
  88     volScalarField r
  89     (
  90         min
  91         (
  92             nuTilda_
  93            /(
  94                max
  95                (
  96                    Stilda,
  97                    dimensionedScalar("SMALL", Stilda.dimensions(), SMALL)
  98                )
  99                *sqr(kappa_*d_)
 100             ),
 101             scalar(10.0)
 102         )
 103     );
 104     r.boundaryField() == 0.0;
 105
 106     volScalarField g(r + Cw2_*(pow6(r) - r));
 107
 108     return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
 109 }
 110
 111
 112 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 113
 114 SpalartAllmaras::SpalartAllmaras
 115 (
 116     const volScalarField& rho,
 117     const volVectorField& U,
 118     const surfaceScalarField& phi,
 119     const basicThermo& thermophysicalModel,
 120     const word& turbulenceModelName,
 121     const word& modelName
 122 )
 123 :
 124     RASModel(modelName, rho, U, phi, thermophysicalModel, turbulenceModelName),
 125
 126     sigmaNut_
 127     (
 128         dimensioned<scalar>::lookupOrAddToDict
 129         (
 130             "sigmaNut",
 131             coeffDict_,
 132             0.66666
 133         )
 134     ),
 135     kappa_
 136     (
 137         dimensioned<scalar>::lookupOrAddToDict
 138         (
 139             "kappa",
 140             coeffDict_,
 141             0.41
 142         )
 143     ),
 144     Prt_
 145     (
 146         dimensioned<scalar>::lookupOrAddToDict
 147         (
 148             "Prt",
 149             coeffDict_,
 150             1.0
 151         )
 152     ),
 153
 154     Cb1_
 155     (
 156         dimensioned<scalar>::lookupOrAddToDict
 157         (
 158             "Cb1",
 159             coeffDict_,
 160             0.1355
 161         )
 162     ),
 163     Cb2_
 164     (
 165         dimensioned<scalar>::lookupOrAddToDict
 166         (
 167             "Cb2",
 168             coeffDict_,
 169             0.622
 170         )
 171     ),
 172     Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
 173     Cw2_
 174     (
 175         dimensioned<scalar>::lookupOrAddToDict
 176         (
 177             "Cw2",
 178             coeffDict_,
 179             0.3
 180         )
 181     ),
 182     Cw3_
 183     (
 184         dimensioned<scalar>::lookupOrAddToDict
 185         (
 186             "Cw3",
 187             coeffDict_,
 188             2.0
 189         )
 190     ),
 191     Cv1_
 192     (
 193         dimensioned<scalar>::lookupOrAddToDict
 194         (
 195             "Cv1",
 196             coeffDict_,
 197             7.1
 198         )
 199     ),
 200     Cv2_
 201     (
 202         dimensioned<scalar>::lookupOrAddToDict
 203         (
 204             "Cv2",
 205             coeffDict_,
 206             5.0
 207         )
 208     ),
 209
 210     nuTilda_
 211     (
 212         IOobject
 213         (
 214             "nuTilda",
 215             runTime_.timeName(),
 216             mesh_,
 217             IOobject::MUST_READ,
 218             IOobject::AUTO_WRITE
 219         ),
 220         mesh_
 221     ),
 222
 223     mut_
 224     (
 225         IOobject
 226         (
 227             "mut",
 228             runTime_.timeName(),
 229             mesh_,
 230             IOobject::MUST_READ,
 231             IOobject::AUTO_WRITE
 232         ),
 233         mesh_
 234     ),
 235
 236     alphat_
 237     (
 238         IOobject
 239         (
 240             "alphat",
 241             runTime_.timeName(),
 242             mesh_,
 243             IOobject::NO_READ,
 244             IOobject::AUTO_WRITE
 245         ),
 246         autoCreateAlphat("alphat", mesh_)
 247     ),
 248
 249     d_(mesh_)
 250 {
 251     alphat_ = mut_/Prt_;
 252     alphat_.correctBoundaryConditions();
 253
 254     printCoeffs();
 255 }
 256
 257
 258 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 259
 260 tmp<volScalarField> SpalartAllmaras::k() const
 261 {
 262     WarningIn("tmp<volScalarField> SpalartAllmaras::k() const")
 263         << "Turbulence kinetic energy not defined for Spalart-Allmaras model. "
 264         << "Returning zero field"
 265         << endl;
 266
 267     return tmp<volScalarField>
 268     (
 269         new volScalarField
 270         (
 271             IOobject
 272             (
 273                 "k",
 274                 runTime_.timeName(),
 275                 mesh_
 276             ),
 277             mesh_,
 278             dimensionedScalar("0", dimensionSet(0, 2, -2, 0, 0), 0)
 279         )
 280     );
 281 }
 282
 283
 284 tmp<volScalarField> SpalartAllmaras::epsilon() const
 285 {
 286     WarningIn("tmp<volScalarField> SpalartAllmaras::epsilon() const")
 287         << "Turbulence kinetic energy dissipation rate not defined for "
 288         << "Spalart-Allmaras model. Returning zero field"
 289         << endl;
 290
 291     return tmp<volScalarField>
 292     (
 293         new volScalarField
 294         (
 295             IOobject
 296             (
 297                 "epslion",
 298                 runTime_.timeName(),
 299                 mesh_
 300             ),
 301             mesh_,
 302             dimensionedScalar("0", dimensionSet(0, 2, -3, 0, 0), 0)
 303         )
 304     );
 305 }
 306
 307
 308 tmp<volSymmTensorField> SpalartAllmaras::R() const
 309 {
 310     return tmp<volSymmTensorField>
 311     (
 312         new volSymmTensorField
 313         (
 314             IOobject
 315             (
 316                 "R",
 317                 runTime_.timeName(),
 318                 mesh_,
 319                 IOobject::NO_READ,
 320                 IOobject::NO_WRITE
 321             ),
 322             ((2.0/3.0)*I)*k() - (mut_/rho_)*dev(twoSymm(fvc::grad(U_)))
 323         )
 324     );
 325 }
 326
 327
 328 tmp<volSymmTensorField> SpalartAllmaras::devRhoReff() const
 329 {
 330     return tmp<volSymmTensorField>
 331     (
 332         new volSymmTensorField
 333         (
 334             IOobject
 335             (
 336                 "devRhoReff",
 337                 runTime_.timeName(),
 338                 mesh_,
 339                 IOobject::NO_READ,
 340                 IOobject::NO_WRITE
 341             ),
 342            -muEff()*dev(twoSymm(fvc::grad(U_)))
 343         )
 344     );
 345 }
 346
 347
 348 tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff(volVectorField& U) const
 349 {
 350     volScalarField muEff_(muEff());
 351
 352     return
 353     (
 354       - fvm::laplacian(muEff_, U)
 355       - fvc::div(muEff_*dev2(T(fvc::grad(U))))
 356     );
 357 }
 358
 359
 360 bool SpalartAllmaras::read()
 361 {
 362     if (RASModel::read())
 363     {
 364         sigmaNut_.readIfPresent(coeffDict());
 365         kappa_.readIfPresent(coeffDict());
 366         Prt_.readIfPresent(coeffDict());
 367
 368         Cb1_.readIfPresent(coeffDict());
 369         Cb2_.readIfPresent(coeffDict());
 370         Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
 371         Cw2_.readIfPresent(coeffDict());
 372         Cw3_.readIfPresent(coeffDict());
 373         Cv1_.readIfPresent(coeffDict());
 374         Cv2_.readIfPresent(coeffDict());
 375
 376         return true;
 377     }
 378     else
 379     {
 380         return false;
 381     }
 382 }
 383
 384
 385 void SpalartAllmaras::correct()
 386 {
 387     if (!turbulence_)
 388     {
 389         // Re-calculate viscosity
 390         mut_ = rho_*nuTilda_*fv1(chi());
 391         mut_.correctBoundaryConditions();
 392
 393         // Re-calculate thermal diffusivity
 394         alphat_ = mut_/Prt_;
 395         alphat_.correctBoundaryConditions();
 396
 397         return;
 398     }
 399
 400     RASModel::correct();
 401
 402     if (mesh_.changing())
 403     {
 404         d_.correct();
 405     }
 406
 407     volScalarField chi(this->chi());
 408     volScalarField fv1(this->fv1(chi));
 409
 410     volScalarField Stilda
 411     (
 412         fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
 413       + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_)
 414     );
 415
 416     tmp<fvScalarMatrix> nuTildaEqn
 417     (
 418         fvm::ddt(rho_, nuTilda_)
 419       + fvm::div(phi_, nuTilda_)
 420       - fvm::laplacian(DnuTildaEff(), nuTilda_)
 421       - Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
 422      ==
 423         Cb1_*rho_*Stilda*nuTilda_
 424       - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_*rho_/sqr(d_), nuTilda_)
 425     );
 426
 427     nuTildaEqn().relax();
 428     solve(nuTildaEqn);
 429     bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
 430     nuTilda_.correctBoundaryConditions();
 431
 432     // Re-calculate viscosity
 433     mut_.internalField() = fv1*nuTilda_.internalField()*rho_.internalField();
 434     mut_.correctBoundaryConditions();
 435
 436     // Re-calculate thermal diffusivity
 437     alphat_ = mut_/Prt_;
 438     alphat_.correctBoundaryConditions();
 439 }
 440
 441
 442 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 443
 444 } // End namespace RASModels
 445 } // End namespace compressible
 446 } // End namespace Foam
 447
 448 // ************************************************************************* //