src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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     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     // Model adjustment, Eric Paterson: boundedness of
  65     // turbulent viscosity.  HJ, 24/Aug/2010
  66     return 1.0 - chi/(1.0 + chi*fv1);
  67     //return 1.0/pow3(scalar(1) + chi/Cv2);
  68 }
  69
  70
  71 tmp<volScalarField> SpalartAllmaras::fv3
  72 (
  73     const volScalarField& chi,
  74     const volScalarField& fv1
  75 ) const
  76 {
  77     volScalarField chiByCv2 = (1/Cv2_)*chi;
  78
  79     return
  80         (scalar(1) + chi*fv1)
  81        *(1/Cv2_)
  82        *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
  83        /pow3(scalar(1) + chiByCv2);
  84 }
  85
  86
  87 tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& Stilda) const
  88 {
  89     volScalarField r = min
  90     (
  91         nuTilda_
  92        /(
  93            max(Stilda, dimensionedScalar("SMALL", Stilda.dimensions(), SMALL))
  94           *sqr(kappa_*d_)
  95         ),
  96         scalar(10.0)
  97     );
  98     r.boundaryField() == 0.0;
  99
 100     volScalarField g = r + Cw2_*(pow6(r) - r);
 101
 102     return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
 103 }
 104
 105
 106 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 107
 108 SpalartAllmaras::SpalartAllmaras
 109 (
 110     const volVectorField& U,
 111     const surfaceScalarField& phi,
 112     transportModel& lamTransportModel
 113 )
 114 :
 115     RASModel(typeName, U, phi, lamTransportModel),
 116
 117     sigmaNut_
 118     (
 119         dimensionedScalar::lookupOrAddToDict
 120         (
 121             "sigmaNut",
 122             coeffDict_,
 123             0.66666
 124         )
 125     ),
 126     kappa_
 127     (
 128         dimensionedScalar::lookupOrAddToDict
 129         (
 130             "kappa",
 131             coeffDict_,
 132             0.41
 133         )
 134     ),
 135
 136     Cb1_
 137     (
 138         dimensionedScalar::lookupOrAddToDict
 139         (
 140             "Cb1",
 141             coeffDict_,
 142             0.1355
 143         )
 144     ),
 145     Cb2_
 146     (
 147         dimensionedScalar::lookupOrAddToDict
 148         (
 149             "Cb2",
 150             coeffDict_,
 151             0.622
 152         )
 153     ),
 154     Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
 155     Cw2_
 156     (
 157         dimensionedScalar::lookupOrAddToDict
 158         (
 159             "Cw2",
 160             coeffDict_,
 161             0.3
 162         )
 163     ),
 164     Cw3_
 165     (
 166         dimensionedScalar::lookupOrAddToDict
 167         (
 168             "Cw3",
 169             coeffDict_,
 170             2.0
 171         )
 172     ),
 173     Cv1_
 174     (
 175         dimensionedScalar::lookupOrAddToDict
 176         (
 177             "Cv1",
 178             coeffDict_,
 179             7.1
 180         )
 181     ),
 182     Cv2_
 183     (
 184         dimensionedScalar::lookupOrAddToDict
 185         (
 186             "Cv2",
 187             coeffDict_,
 188             5.0
 189         )
 190     ),
 191
 192     nuTilda_
 193     (
 194         IOobject
 195         (
 196             "nuTilda",
 197             runTime_.timeName(),
 198             U_.db(),
 199             IOobject::MUST_READ,
 200             IOobject::AUTO_WRITE
 201         ),
 202         mesh_
 203     ),
 204
 205     nut_
 206     (
 207         IOobject
 208         (
 209             "nut",
 210             runTime_.timeName(),
 211             U_.db(),
 212             IOobject::MUST_READ,
 213             IOobject::AUTO_WRITE
 214         ),
 215         mesh_
 216     ),
 217
 218     d_(mesh_)
 219 {
 220     printCoeffs();
 221 }
 222
 223
 224 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 225
 226 tmp<volScalarField> SpalartAllmaras::DnuTildaEff() const
 227 {
 228     return tmp<volScalarField>
 229     (
 230         new volScalarField("DnuTildaEff", nuTilda_/sigmaNut_ + nu())
 231     );
 232 }
 233
 234
 235 tmp<volScalarField> SpalartAllmaras::k() const
 236 {
 237     return tmp<volScalarField>
 238     (
 239         new volScalarField
 240         (
 241             IOobject
 242             (
 243                 "k",
 244                 runTime_.timeName(),
 245                 U_.db()
 246             ),
 247             mesh_,
 248             dimensionedScalar("0", dimensionSet(0, 2, -2, 0, 0), 0)
 249         )
 250     );
 251 }
 252
 253
 254 tmp<volScalarField> SpalartAllmaras::epsilon() const
 255 {
 256     return tmp<volScalarField>
 257     (
 258         new volScalarField
 259         (
 260             IOobject
 261             (
 262                 "epsilon",
 263                 runTime_.timeName(),
 264                 U_.db()
 265             ),
 266             mesh_,
 267             dimensionedScalar("0", dimensionSet(0, 2, -3, 0, 0), 0)
 268         )
 269     );
 270 }
 271
 272
 273 tmp<volSymmTensorField> SpalartAllmaras::R() const
 274 {
 275     return tmp<volSymmTensorField>
 276     (
 277         new volSymmTensorField
 278         (
 279             IOobject
 280             (
 281                 "R",
 282                 runTime_.timeName(),
 283                 U_.db(),
 284                 IOobject::NO_READ,
 285                 IOobject::NO_WRITE
 286             ),
 287             ((2.0/3.0)*I)*k() - nut()*twoSymm(fvc::grad(U_))
 288         )
 289     );
 290 }
 291
 292
 293 tmp<volSymmTensorField> SpalartAllmaras::devReff() const
 294 {
 295     return tmp<volSymmTensorField>
 296     (
 297         new volSymmTensorField
 298         (
 299             IOobject
 300             (
 301                 "devRhoReff",
 302                 runTime_.timeName(),
 303                 U_.db(),
 304                 IOobject::NO_READ,
 305                 IOobject::NO_WRITE
 306             ),
 307            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 308         )
 309     );
 310 }
 311
 312
 313 tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
 314 {
 315     volScalarField nuEff_ = nuEff();
 316
 317     return
 318     (
 319       - fvm::laplacian(nuEff_, U)
 320       - fvc::div(nuEff_*dev(fvc::grad(U)().T()))
 321     );
 322 }
 323
 324
 325 bool SpalartAllmaras::read()
 326 {
 327     if (RASModel::read())
 328     {
 329         sigmaNut_.readIfPresent(coeffDict());
 330
 331         Cb1_.readIfPresent(coeffDict());
 332         Cb2_.readIfPresent(coeffDict());
 333         Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
 334         Cw2_.readIfPresent(coeffDict());
 335         Cw3_.readIfPresent(coeffDict());
 336         Cv1_.readIfPresent(coeffDict());
 337         Cv2_.readIfPresent(coeffDict());
 338
 339         return true;
 340     }
 341     else
 342     {
 343         return false;
 344     }
 345 }
 346
 347
 348 void SpalartAllmaras::correct()
 349 {
 350     if (mesh_.changing())
 351     {
 352         bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
 353     }
 354
 355     RASModel::correct();
 356
 357     if (!turbulence_)
 358     {
 359         return;
 360     }
 361
 362     if (mesh_.changing())
 363     {
 364         d_.correct();
 365     }
 366
 367     volScalarField chi = this->chi();
 368     volScalarField fv1 = this->fv1(chi);
 369
 370     volScalarField Stilda =
 371 //        fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
 372         sqrt(2.0)*mag(skew(fvc::grad(U_)))
 373       + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_);
 374
 375     tmp<fvScalarMatrix> nuTildaEqn
 376     (
 377         fvm::ddt(nuTilda_)
 378       + fvm::div(phi_, nuTilda_)
 379       + fvm::SuSp(-fvc::div(phi_), nuTilda_)
 380       - fvm::laplacian(DnuTildaEff(), nuTilda_)
 381       - Cb2_/sigmaNut_*magSqr(fvc::grad(nuTilda_))
 382      ==
 383         Cb1_*Stilda*nuTilda_
 384       - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_/sqr(d_), nuTilda_)
 385     );
 386
 387     nuTildaEqn().relax();
 388     solve(nuTildaEqn);
 389     bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
 390     nuTilda_.correctBoundaryConditions();
 391
 392     nut_.internalField() = fv1*nuTilda_.internalField();
 393     nut_ = min(nut_, nuRatio()*nu());
 394     nut_.correctBoundaryConditions();
 395 }
 396
 397
 398 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 399
 400 } // End namespace RASModels
 401 } // End namespace incompressible
 402 } // End namespace Foam
 403
 404 // ************************************************************************* //