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