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