src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 #include "backwardsCompatibilityWallFunctions.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36 namespace compressible
  37 {
  38 namespace RASModels
  39 {
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 defineTypeNameAndDebug(SpalartAllmaras, 0);
  44 addToRunTimeSelectionTable(RASModel, SpalartAllmaras, dictionary);
  45
  46 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  47
  48 tmp<volScalarField> SpalartAllmaras::chi() const
  49 {
  50     return rho_*nuTilda_/mu();
  51 }
  52
  53
  54 tmp<volScalarField> SpalartAllmaras::fv1(const volScalarField& chi) const
  55 {
  56     volScalarField chi3 = pow3(chi);
  57     return chi3/(chi3 + pow3(Cv1_));
  58 }
  59
  60
  61 tmp<volScalarField> SpalartAllmaras::fv2
  62 (
  63     const volScalarField& chi,
  64     const volScalarField& fv1
  65 ) const
  66 {
  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 volScalarField& rho,
 111     const volVectorField& U,
 112     const surfaceScalarField& phi,
 113     const basicThermo& thermophysicalModel
 114 )
 115 :
 116     RASModel(typeName, rho, U, phi, thermophysicalModel),
 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     Prt_
 137     (
 138         dimensioned<scalar>::lookupOrAddToDict
 139         (
 140             "Prt",
 141             coeffDict_,
 142             1.0
 143         )
 144     ),
 145
 146     Cb1_
 147     (
 148         dimensioned<scalar>::lookupOrAddToDict
 149         (
 150             "Cb1",
 151             coeffDict_,
 152             0.1355
 153         )
 154     ),
 155     Cb2_
 156     (
 157         dimensioned<scalar>::lookupOrAddToDict
 158         (
 159             "Cb2",
 160             coeffDict_,
 161             0.622
 162         )
 163     ),
 164     Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
 165     Cw2_
 166     (
 167         dimensioned<scalar>::lookupOrAddToDict
 168         (
 169             "Cw2",
 170             coeffDict_,
 171             0.3
 172         )
 173     ),
 174     Cw3_
 175     (
 176         dimensioned<scalar>::lookupOrAddToDict
 177         (
 178             "Cw3",
 179             coeffDict_,
 180             2.0
 181         )
 182     ),
 183     Cv1_
 184     (
 185         dimensioned<scalar>::lookupOrAddToDict
 186         (
 187             "Cv1",
 188             coeffDict_,
 189             7.1
 190         )
 191     ),
 192     Cv2_
 193     (
 194         dimensioned<scalar>::lookupOrAddToDict
 195         (
 196             "Cv2",
 197             coeffDict_,
 198             5.0
 199         )
 200     ),
 201
 202     nuTilda_
 203     (
 204         IOobject
 205         (
 206             "nuTilda",
 207             runTime_.timeName(),
 208             mesh_,
 209             IOobject::MUST_READ,
 210             IOobject::AUTO_WRITE
 211         ),
 212         mesh_
 213     ),
 214
 215     mut_
 216     (
 217         IOobject
 218         (
 219             "mut",
 220             runTime_.timeName(),
 221             mesh_,
 222             IOobject::MUST_READ,
 223             IOobject::AUTO_WRITE
 224         ),
 225         mesh_
 226     ),
 227
 228     alphat_
 229     (
 230         IOobject
 231         (
 232             "alphat",
 233             runTime_.timeName(),
 234             mesh_,
 235             IOobject::NO_READ,
 236             IOobject::AUTO_WRITE
 237         ),
 238         autoCreateAlphat("alphat", mesh_)
 239     ),
 240
 241     d_(mesh_)
 242 {
 243     alphat_ = mut_/Prt_;
 244     alphat_.correctBoundaryConditions();
 245
 246     printCoeffs();
 247 }
 248
 249
 250 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 251
 252 tmp<volSymmTensorField> SpalartAllmaras::R() const
 253 {
 254     return tmp<volSymmTensorField>
 255     (
 256         new volSymmTensorField
 257         (
 258             IOobject
 259             (
 260                 "R",
 261                 runTime_.timeName(),
 262                 mesh_,
 263                 IOobject::NO_READ,
 264                 IOobject::NO_WRITE
 265             ),
 266             ((2.0/3.0)*I)*k() - (mut_/rho_)*dev(twoSymm(fvc::grad(U_)))
 267         )
 268     );
 269 }
 270
 271
 272 tmp<volSymmTensorField> SpalartAllmaras::devRhoReff() const
 273 {
 274     return tmp<volSymmTensorField>
 275     (
 276         new volSymmTensorField
 277         (
 278             IOobject
 279             (
 280                 "devRhoReff",
 281                 runTime_.timeName(),
 282                 mesh_,
 283                 IOobject::NO_READ,
 284                 IOobject::NO_WRITE
 285             ),
 286            -muEff()*dev(twoSymm(fvc::grad(U_)))
 287         )
 288     );
 289 }
 290
 291
 292 tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff(volVectorField& U) const
 293 {
 294     volScalarField muEff_ = muEff();
 295
 296     return
 297     (
 298       - fvm::laplacian(muEff_, U)
 299       - fvc::div(muEff_*dev2(fvc::grad(U)().T()))
 300     );
 301 }
 302
 303
 304 bool SpalartAllmaras::read()
 305 {
 306     if (RASModel::read())
 307     {
 308         sigmaNut_.readIfPresent(coeffDict());
 309         kappa_.readIfPresent(coeffDict());
 310         Prt_.readIfPresent(coeffDict());
 311
 312         Cb1_.readIfPresent(coeffDict());
 313         Cb2_.readIfPresent(coeffDict());
 314         Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
 315         Cw2_.readIfPresent(coeffDict());
 316         Cw3_.readIfPresent(coeffDict());
 317         Cv1_.readIfPresent(coeffDict());
 318         Cv2_.readIfPresent(coeffDict());
 319
 320         return true;
 321     }
 322     else
 323     {
 324         return false;
 325     }
 326 }
 327
 328
 329 void SpalartAllmaras::correct()
 330 {
 331     if (!turbulence_)
 332     {
 333         // Re-calculate viscosity
 334         mut_ = rho_*nuTilda_*fv1(chi());
 335         mut_.correctBoundaryConditions();
 336
 337         // Re-calculate thermal diffusivity
 338         alphat_ = mut_/Prt_;
 339         alphat_.correctBoundaryConditions();
 340
 341         return;
 342     }
 343
 344     RASModel::correct();
 345
 346     if (mesh_.changing())
 347     {
 348         d_.correct();
 349     }
 350
 351     volScalarField chi = this->chi();
 352     volScalarField fv1 = this->fv1(chi);
 353
 354     volScalarField Stilda =
 355         fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
 356       + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_);
 357
 358     tmp<fvScalarMatrix> nuTildaEqn
 359     (
 360         fvm::ddt(rho_, nuTilda_)
 361       + fvm::div(phi_, nuTilda_)
 362       - fvm::laplacian(DnuTildaEff(), nuTilda_)
 363       - Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
 364      ==
 365         Cb1_*rho_*Stilda*nuTilda_
 366       - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_*rho_/sqr(d_), nuTilda_)
 367     );
 368
 369     nuTildaEqn().relax();
 370     solve(nuTildaEqn);
 371     bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
 372     nuTilda_.correctBoundaryConditions();
 373
 374     // Re-calculate viscosity
 375     mut_.internalField() = fv1*nuTilda_.internalField()*rho_.internalField();
 376     mut_.correctBoundaryConditions();
 377
 378     // Re-calculate thermal diffusivity
 379     alphat_ = mut_/Prt_;
 380     alphat_.correctBoundaryConditions();
 381 }
 382
 383
 384 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 385
 386 } // End namespace RASModels
 387 } // End namespace compressible
 388 } // End namespace Foam
 389
 390 // ************************************************************************* //