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