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