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