src/turbulenceModels/incompressible/RAS/kOmega/kOmega.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 "kOmega.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(kOmega, 0);
  43 addToRunTimeSelectionTable(RASModel, kOmega, dictionary);
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 kOmega::kOmega
  48 (
  49     const volVectorField& U,
  50     const surfaceScalarField& phi,
  51     transportModel& transport,
  52     const word& turbulenceModelName,
  53     const word& modelName
  54 )
  55 :
  56     RASModel(modelName, U, phi, transport, turbulenceModelName),
  57
  58     Cmu_
  59     (
  60         dimensioned<scalar>::lookupOrAddToDict
  61         (
  62             "betaStar",
  63             coeffDict_,
  64             0.09
  65         )
  66     ),
  67     beta_
  68     (
  69         dimensioned<scalar>::lookupOrAddToDict
  70         (
  71             "beta",
  72             coeffDict_,
  73             0.072
  74         )
  75     ),
  76     alpha_
  77     (
  78         dimensioned<scalar>::lookupOrAddToDict
  79         (
  80             "alpha",
  81             coeffDict_,
  82             0.52
  83         )
  84     ),
  85     alphaK_
  86     (
  87         dimensioned<scalar>::lookupOrAddToDict
  88         (
  89             "alphaK",
  90             coeffDict_,
  91             0.5
  92         )
  93     ),
  94     alphaOmega_
  95     (
  96         dimensioned<scalar>::lookupOrAddToDict
  97         (
  98             "alphaOmega",
  99             coeffDict_,
 100             0.5
 101         )
 102     ),
 103
 104     k_
 105     (
 106         IOobject
 107         (
 108             "k",
 109             runTime_.timeName(),
 110             mesh_,
 111             IOobject::NO_READ,
 112             IOobject::AUTO_WRITE
 113         ),
 114         autoCreateK("k", mesh_)
 115     ),
 116     omega_
 117     (
 118         IOobject
 119         (
 120             "omega",
 121             runTime_.timeName(),
 122             mesh_,
 123             IOobject::NO_READ,
 124             IOobject::AUTO_WRITE
 125         ),
 126         autoCreateOmega("omega", mesh_)
 127     ),
 128     nut_
 129     (
 130         IOobject
 131         (
 132             "nut",
 133             runTime_.timeName(),
 134             mesh_,
 135             IOobject::NO_READ,
 136             IOobject::AUTO_WRITE
 137         ),
 138         autoCreateNut("nut", mesh_)
 139     )
 140 {
 141     bound(k_, kMin_);
 142     bound(omega_, omegaMin_);
 143
 144     nut_ = k_/omega_;
 145     nut_.correctBoundaryConditions();
 146
 147     printCoeffs();
 148 }
 149
 150
 151 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 152
 153 tmp<volSymmTensorField> kOmega::R() const
 154 {
 155     return tmp<volSymmTensorField>
 156     (
 157         new volSymmTensorField
 158         (
 159             IOobject
 160             (
 161                 "R",
 162                 runTime_.timeName(),
 163                 mesh_,
 164                 IOobject::NO_READ,
 165                 IOobject::NO_WRITE
 166             ),
 167             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 168             k_.boundaryField().types()
 169         )
 170     );
 171 }
 172
 173
 174 tmp<volSymmTensorField> kOmega::devReff() const
 175 {
 176     return tmp<volSymmTensorField>
 177     (
 178         new volSymmTensorField
 179         (
 180             IOobject
 181             (
 182                 "devRhoReff",
 183                 runTime_.timeName(),
 184                 mesh_,
 185                 IOobject::NO_READ,
 186                 IOobject::NO_WRITE
 187             ),
 188            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 189         )
 190     );
 191 }
 192
 193
 194 tmp<fvVectorMatrix> kOmega::divDevReff(volVectorField& U) const
 195 {
 196     return
 197     (
 198       - fvm::laplacian(nuEff(), U)
 199       - fvc::div(nuEff()*dev(T(fvc::grad(U))))
 200     );
 201 }
 202
 203
 204 bool kOmega::read()
 205 {
 206     if (RASModel::read())
 207     {
 208         Cmu_.readIfPresent(coeffDict());
 209         beta_.readIfPresent(coeffDict());
 210         alphaK_.readIfPresent(coeffDict());
 211         alphaOmega_.readIfPresent(coeffDict());
 212
 213         return true;
 214     }
 215     else
 216     {
 217         return false;
 218     }
 219 }
 220
 221
 222 void kOmega::correct()
 223 {
 224     RASModel::correct();
 225
 226     if (!turbulence_)
 227     {
 228         return;
 229     }
 230
 231     volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_))));
 232
 233     // Update omega and G at the wall
 234     omega_.boundaryField().updateCoeffs();
 235
 236     // Turbulence specific dissipation rate equation
 237     tmp<fvScalarMatrix> omegaEqn
 238     (
 239         fvm::ddt(omega_)
 240       + fvm::div(phi_, omega_)
 241       - fvm::Sp(fvc::div(phi_), omega_)
 242       - fvm::laplacian(DomegaEff(), omega_)
 243      ==
 244         alpha_*G*omega_/k_
 245       - fvm::Sp(beta_*omega_, omega_)
 246     );
 247
 248     omegaEqn().relax();
 249
 250     omegaEqn().boundaryManipulate(omega_.boundaryField());
 251
 252     solve(omegaEqn);
 253     bound(omega_, omegaMin_);
 254
 255
 256     // Turbulent kinetic energy equation
 257     tmp<fvScalarMatrix> kEqn
 258     (
 259         fvm::ddt(k_)
 260       + fvm::div(phi_, k_)
 261       - fvm::Sp(fvc::div(phi_), k_)
 262       - fvm::laplacian(DkEff(), k_)
 263      ==
 264         G
 265       - fvm::Sp(Cmu_*omega_, k_)
 266     );
 267
 268     kEqn().relax();
 269     solve(kEqn);
 270     bound(k_, kMin_);
 271
 272
 273     // Re-calculate viscosity
 274     nut_ = k_/omega_;
 275     nut_.correctBoundaryConditions();
 276 }
 277
 278
 279 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 280
 281 } // End namespace RASModels
 282 } // End namespace incompressible
 283 } // End namespace Foam
 284
 285 // ************************************************************************* //