src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  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 "qZeta.H"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 #include "backwardsCompatibilityWallFunctions.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35 namespace incompressible
  36 {
  37 namespace RASModels
  38 {
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(qZeta, 0);
  43 addToRunTimeSelectionTable(RASModel, qZeta, dictionary);
  44
  45 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  46
  47 tmp<volScalarField> qZeta::fMu() const
  48 {
  49     const volScalarField Rt(q_*k_/(2.0*nu()*zeta_));
  50
  51     if (anisotropic_)
  52     {
  53         return exp((-scalar(2.5) + Rt/20.0)/pow(scalar(1) + Rt/130.0, 3.0));
  54     }
  55     else
  56     {
  57         return
  58             exp(-6.0/sqr(scalar(1) + Rt/50.0))
  59            *(scalar(1) + 3.0*exp(-Rt/10.0));
  60     }
  61 }
  62
  63
  64 tmp<volScalarField> qZeta::f2() const
  65 {
  66     tmp<volScalarField> Rt = q_*k_/(2.0*nu()*zeta_);
  67     return scalar(1) - 0.3*exp(-sqr(Rt));
  68 }
  69
  70
  71 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  72
  73 qZeta::qZeta
  74 (
  75     const volVectorField& U,
  76     const surfaceScalarField& phi,
  77     transportModel& transport,
  78     const word& turbulenceModelName,
  79     const word& modelName
  80 )
  81 :
  82     RASModel(modelName, U, phi, transport, turbulenceModelName),
  83
  84     Cmu_
  85     (
  86         dimensioned<scalar>::lookupOrAddToDict
  87         (
  88             "Cmu",
  89             coeffDict_,
  90             0.09
  91         )
  92     ),
  93     C1_
  94     (
  95         dimensioned<scalar>::lookupOrAddToDict
  96         (
  97             "C1",
  98             coeffDict_,
  99             1.44
 100         )
 101     ),
 102     C2_
 103     (
 104         dimensioned<scalar>::lookupOrAddToDict
 105         (
 106             "C2",
 107             coeffDict_,
 108             1.92
 109         )
 110     ),
 111     sigmaZeta_
 112     (
 113         dimensioned<scalar>::lookupOrAddToDict
 114         (
 115             "sigmaZeta",
 116             coeffDict_,
 117             1.3
 118         )
 119     ),
 120     anisotropic_
 121     (
 122         Switch::lookupOrAddToDict
 123         (
 124             "anisotropic",
 125             coeffDict_,
 126             false
 127         )
 128     ),
 129
 130     qMin_("qMin", dimVelocity, SMALL),
 131     zetaMin_("zetaMin", dimVelocity/dimTime, SMALL),
 132
 133     k_
 134     (
 135         IOobject
 136         (
 137             "k",
 138             runTime_.timeName(),
 139             mesh_,
 140             IOobject::MUST_READ,
 141             IOobject::AUTO_WRITE
 142         ),
 143         mesh_
 144     ),
 145
 146     epsilon_
 147     (
 148         IOobject
 149         (
 150             "epsilon",
 151             runTime_.timeName(),
 152             mesh_,
 153             IOobject::MUST_READ,
 154             IOobject::AUTO_WRITE
 155         ),
 156         mesh_
 157     ),
 158
 159     q_
 160     (
 161         IOobject
 162         (
 163             "q",
 164             runTime_.timeName(),
 165             mesh_,
 166             IOobject::NO_READ,
 167             IOobject::AUTO_WRITE
 168         ),
 169         sqrt(k_),
 170         k_.boundaryField().types()
 171     ),
 172
 173     zeta_
 174     (
 175         IOobject
 176         (
 177             "zeta",
 178             runTime_.timeName(),
 179             mesh_,
 180             IOobject::NO_READ,
 181             IOobject::AUTO_WRITE
 182         ),
 183         epsilon_/(2.0*bound(q_, qMin_)),
 184         epsilon_.boundaryField().types()
 185     ),
 186
 187     nut_
 188     (
 189         IOobject
 190         (
 191             "nut",
 192             runTime_.timeName(),
 193             mesh_,
 194             IOobject::NO_READ,
 195             IOobject::AUTO_WRITE
 196         ),
 197         autoCreateNut("nut", mesh_)
 198     )
 199 {
 200     bound(k_, kMin_);
 201     bound(epsilon_, epsilonMin_);
 202     // already bounded: bound(q_, qMin_);
 203     bound(zeta_, zetaMin_);
 204
 205     nut_ = Cmu_*fMu()*sqr(k_)/epsilon_;
 206     nut_.correctBoundaryConditions();
 207
 208     printCoeffs();
 209 }
 210
 211
 212 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 213
 214 tmp<volSymmTensorField> qZeta::R() const
 215 {
 216     return tmp<volSymmTensorField>
 217     (
 218         new volSymmTensorField
 219         (
 220             IOobject
 221             (
 222                 "R",
 223                 runTime_.timeName(),
 224                 mesh_,
 225                 IOobject::NO_READ,
 226                 IOobject::NO_WRITE
 227             ),
 228             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 229             k_.boundaryField().types()
 230         )
 231     );
 232 }
 233
 234
 235 tmp<volSymmTensorField> qZeta::devReff() const
 236 {
 237     return tmp<volSymmTensorField>
 238     (
 239         new volSymmTensorField
 240         (
 241             IOobject
 242             (
 243                 "devRhoReff",
 244                 runTime_.timeName(),
 245                 mesh_,
 246                 IOobject::NO_READ,
 247                 IOobject::NO_WRITE
 248             ),
 249            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 250         )
 251     );
 252 }
 253
 254
 255 tmp<fvVectorMatrix> qZeta::divDevReff(volVectorField& U) const
 256 {
 257     return
 258     (
 259       - fvm::laplacian(nuEff(), U)
 260       - fvc::div(nuEff()*dev(T(fvc::grad(U))))
 261     );
 262 }
 263
 264
 265 bool qZeta::read()
 266 {
 267     if (RASModel::read())
 268     {
 269         Cmu_.readIfPresent(coeffDict());
 270         C1_.readIfPresent(coeffDict());
 271         C2_.readIfPresent(coeffDict());
 272         sigmaZeta_.readIfPresent(coeffDict());
 273         anisotropic_.readIfPresent("anisotropic", coeffDict());
 274
 275         qMin_.readIfPresent(*this);
 276         zetaMin_.readIfPresent(*this);
 277
 278         return true;
 279     }
 280     else
 281     {
 282         return false;
 283     }
 284 }
 285
 286
 287 void qZeta::correct()
 288 {
 289     RASModel::correct();
 290
 291     if (!turbulence_)
 292     {
 293         return;
 294     }
 295
 296     tmp<volScalarField> S2 = 2*magSqr(symm(fvc::grad(U_)));
 297
 298     volScalarField G("RASModel::G", nut_/(2.0*q_)*S2);
 299     const volScalarField E(nu()*nut_/q_*fvc::magSqrGradGrad(U_));
 300
 301
 302     // Zeta equation
 303
 304     tmp<fvScalarMatrix> zetaEqn
 305     (
 306         fvm::ddt(zeta_)
 307       + fvm::div(phi_, zeta_)
 308       - fvm::laplacian(DzetaEff(), zeta_)
 309      ==
 310         (2.0*C1_ - 1)*G*zeta_/q_
 311       - fvm::Sp((2.0*C2_ - dimensionedScalar(1.0))*f2()*zeta_/q_, zeta_)
 312       + E
 313     );
 314
 315     zetaEqn().relax();
 316     solve(zetaEqn);
 317     bound(zeta_, zetaMin_);
 318
 319
 320     // q equation
 321
 322     tmp<fvScalarMatrix> qEqn
 323     (
 324         fvm::ddt(q_)
 325       + fvm::div(phi_, q_)
 326       - fvm::laplacian(DqEff(), q_)
 327      ==
 328         G - fvm::Sp(zeta_/q_, q_)
 329     );
 330
 331     qEqn().relax();
 332     solve(qEqn);
 333     bound(q_, qMin_);
 334
 335
 336     // Re-calculate k and epsilon
 337     k_ = sqr(q_);
 338     k_.correctBoundaryConditions();
 339
 340     epsilon_ = 2*q_*zeta_;
 341     epsilon_.correctBoundaryConditions();
 342
 343
 344     // Re-calculate viscosity
 345     nut_ = Cmu_*fMu()*sqr(k_)/epsilon_;
 346     nut_.correctBoundaryConditions();
 347 }
 348
 349
 350 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 351
 352 } // End namespace RASModels
 353 } // End namespace incompressible
 354 } // End namespace Foam
 355
 356 // ************************************************************************* //