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

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