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) 1991-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 = min
  89     (
  90         nuTilda_
  91        /(
  92            max(Stilda, dimensionedScalar("SMALL", Stilda.dimensions(), SMALL))
  93            *sqr(kappa_*d_)
  94         ),
  95         scalar(10.0)
  96     );
  97     r.boundaryField() == 0.0;
  98
  99     volScalarField g = r + Cw2_*(pow6(r) - r);
 100
 101     return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
 102 }
 103
 104
 105 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 106
 107 SpalartAllmaras::SpalartAllmaras
 108 (
 109     const volScalarField& rho,
 110     const volVectorField& U,
 111     const surfaceScalarField& phi,
 112     const basicThermo& thermophysicalModel
 113 )
 114 :
 115     RASModel(typeName, rho, U, phi, thermophysicalModel),
 116
 117     sigmaNut_
 118     (
 119         dimensioned<scalar>::lookupOrAddToDict
 120         (
 121             "sigmaNut",
 122             coeffDict_,
 123             0.66666
 124         )
 125     ),
 126     kappa_
 127     (
 128         dimensioned<scalar>::lookupOrAddToDict
 129         (
 130             "kappa",
 131             coeffDict_,
 132             0.41
 133         )
 134     ),
 135     Prt_
 136     (
 137         dimensioned<scalar>::lookupOrAddToDict
 138         (
 139             "Prt",
 140             coeffDict_,
 141             1.0
 142         )
 143     ),
 144
 145     Cb1_
 146     (
 147         dimensioned<scalar>::lookupOrAddToDict
 148         (
 149             "Cb1",
 150             coeffDict_,
 151             0.1355
 152         )
 153     ),
 154     Cb2_
 155     (
 156         dimensioned<scalar>::lookupOrAddToDict
 157         (
 158             "Cb2",
 159             coeffDict_,
 160             0.622
 161         )
 162     ),
 163     Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
 164     Cw2_
 165     (
 166         dimensioned<scalar>::lookupOrAddToDict
 167         (
 168             "Cw2",
 169             coeffDict_,
 170             0.3
 171         )
 172     ),
 173     Cw3_
 174     (
 175         dimensioned<scalar>::lookupOrAddToDict
 176         (
 177             "Cw3",
 178             coeffDict_,
 179             2.0
 180         )
 181     ),
 182     Cv1_
 183     (
 184         dimensioned<scalar>::lookupOrAddToDict
 185         (
 186             "Cv1",
 187             coeffDict_,
 188             7.1
 189         )
 190     ),
 191     Cv2_
 192     (
 193         dimensioned<scalar>::lookupOrAddToDict
 194         (
 195             "Cv2",
 196             coeffDict_,
 197             5.0
 198         )
 199     ),
 200
 201     nuTilda_
 202     (
 203         IOobject
 204         (
 205             "nuTilda",
 206             runTime_.timeName(),
 207             mesh_,
 208             IOobject::MUST_READ,
 209             IOobject::AUTO_WRITE
 210         ),
 211         mesh_
 212     ),
 213
 214     mut_
 215     (
 216         IOobject
 217         (
 218             "mut",
 219             runTime_.timeName(),
 220             mesh_,
 221             IOobject::MUST_READ,
 222             IOobject::AUTO_WRITE
 223         ),
 224         mesh_
 225     ),
 226
 227     alphat_
 228     (
 229         IOobject
 230         (
 231             "alphat",
 232             runTime_.timeName(),
 233             mesh_,
 234             IOobject::NO_READ,
 235             IOobject::AUTO_WRITE
 236         ),
 237         autoCreateAlphat("alphat", mesh_)
 238     ),
 239
 240     d_(mesh_)
 241 {
 242     alphat_ = mut_/Prt_;
 243     alphat_.correctBoundaryConditions();
 244
 245     printCoeffs();
 246 }
 247
 248
 249 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 250
 251 tmp<volSymmTensorField> SpalartAllmaras::R() const
 252 {
 253     return tmp<volSymmTensorField>
 254     (
 255         new volSymmTensorField
 256         (
 257             IOobject
 258             (
 259                 "R",
 260                 runTime_.timeName(),
 261                 mesh_,
 262                 IOobject::NO_READ,
 263                 IOobject::NO_WRITE
 264             ),
 265             ((2.0/3.0)*I)*k() - (mut_/rho_)*dev(twoSymm(fvc::grad(U_)))
 266         )
 267     );
 268 }
 269
 270
 271 tmp<volSymmTensorField> SpalartAllmaras::devRhoReff() const
 272 {
 273     return tmp<volSymmTensorField>
 274     (
 275         new volSymmTensorField
 276         (
 277             IOobject
 278             (
 279                 "devRhoReff",
 280                 runTime_.timeName(),
 281                 mesh_,
 282                 IOobject::NO_READ,
 283                 IOobject::NO_WRITE
 284             ),
 285            -muEff()*dev(twoSymm(fvc::grad(U_)))
 286         )
 287     );
 288 }
 289
 290
 291 tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff(volVectorField& U) const
 292 {
 293     volScalarField muEff_ = muEff();
 294
 295     return
 296     (
 297       - fvm::laplacian(muEff_, U)
 298       - fvc::div(muEff_*dev2(fvc::grad(U)().T()))
 299     );
 300 }
 301
 302
 303 bool SpalartAllmaras::read()
 304 {
 305     if (RASModel::read())
 306     {
 307         sigmaNut_.readIfPresent(coeffDict());
 308         kappa_.readIfPresent(coeffDict());
 309         Prt_.readIfPresent(coeffDict());
 310
 311         Cb1_.readIfPresent(coeffDict());
 312         Cb2_.readIfPresent(coeffDict());
 313         Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
 314         Cw2_.readIfPresent(coeffDict());
 315         Cw3_.readIfPresent(coeffDict());
 316         Cv1_.readIfPresent(coeffDict());
 317         Cv2_.readIfPresent(coeffDict());
 318
 319         return true;
 320     }
 321     else
 322     {
 323         return false;
 324     }
 325 }
 326
 327
 328 void SpalartAllmaras::correct()
 329 {
 330     if (!turbulence_)
 331     {
 332         // Re-calculate viscosity
 333         mut_ = rho_*nuTilda_*fv1(chi());
 334         mut_.correctBoundaryConditions();
 335
 336         // Re-calculate thermal diffusivity
 337         alphat_ = mut_/Prt_;
 338         alphat_.correctBoundaryConditions();
 339
 340         return;
 341     }
 342
 343     RASModel::correct();
 344
 345     if (mesh_.changing())
 346     {
 347         d_.correct();
 348     }
 349
 350     volScalarField chi = this->chi();
 351     volScalarField fv1 = this->fv1(chi);
 352
 353     volScalarField Stilda =
 354         fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
 355       + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_);
 356
 357     tmp<fvScalarMatrix> nuTildaEqn
 358     (
 359         fvm::ddt(rho_, nuTilda_)
 360       + fvm::div(phi_, nuTilda_)
 361       - fvm::laplacian(DnuTildaEff(), nuTilda_)
 362       - Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
 363      ==
 364         Cb1_*rho_*Stilda*nuTilda_
 365       - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_*rho_/sqr(d_), nuTilda_)
 366     );
 367
 368     nuTildaEqn().relax();
 369     solve(nuTildaEqn);
 370     bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
 371     nuTilda_.correctBoundaryConditions();
 372
 373     // Re-calculate viscosity
 374     mut_.internalField() = fv1*nuTilda_.internalField()*rho_.internalField();
 375     mut_.correctBoundaryConditions();
 376
 377     // Re-calculate thermal diffusivity
 378     alphat_ = mut_/Prt_;
 379     alphat_.correctBoundaryConditions();
 380 }
 381
 382
 383 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 384
 385 } // End namespace RASModels
 386 } // End namespace compressible
 387 } // End namespace Foam
 388
 389 // ************************************************************************* //