src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.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 "LamBremhorstKE.H"
  28 #include "addToRunTimeSelectionTable.H"
  29
  30 #include "backwardsCompatibilityWallFunctions.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36 namespace incompressible
  37 {
  38 namespace RASModels
  39 {
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 defineTypeNameAndDebug(LamBremhorstKE, 0);
  44 addToRunTimeSelectionTable(RASModel, LamBremhorstKE, dictionary);
  45
  46 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  47
  48 LamBremhorstKE::LamBremhorstKE
  49 (
  50     const volVectorField& U,
  51     const surfaceScalarField& phi,
  52     transportModel& lamTransportModel
  53 )
  54 :
  55     RASModel(typeName, U, phi, lamTransportModel),
  56
  57     Cmu_
  58     (
  59         dimensioned<scalar>::lookupOrAddToDict
  60         (
  61             "Cmu",
  62             coeffDict_,
  63             0.09
  64         )
  65     ),
  66     C1_
  67     (
  68         dimensioned<scalar>::lookupOrAddToDict
  69         (
  70             "C1",
  71             coeffDict_,
  72             1.44
  73         )
  74     ),
  75     C2_
  76     (
  77         dimensioned<scalar>::lookupOrAddToDict
  78         (
  79             "C2",
  80             coeffDict_,
  81             1.92
  82         )
  83     ),
  84     sigmaEps_
  85     (
  86         dimensioned<scalar>::lookupOrAddToDict
  87         (
  88             "alphaEps",
  89             coeffDict_,
  90             1.3
  91         )
  92     ),
  93
  94     k_
  95     (
  96         IOobject
  97         (
  98             "k",
  99             runTime_.timeName(),
 100             mesh_,
 101             IOobject::MUST_READ,
 102             IOobject::AUTO_WRITE
 103         ),
 104         mesh_
 105     ),
 106
 107     epsilon_
 108     (
 109         IOobject
 110         (
 111             "epsilon",
 112             runTime_.timeName(),
 113             mesh_,
 114             IOobject::MUST_READ,
 115             IOobject::AUTO_WRITE
 116         ),
 117         autoCreateEpsilon("epsilon", mesh_)
 118     ),
 119
 120     y_(mesh_),
 121
 122     Rt_(sqr(k_)/(nu()*epsilon_)),
 123
 124     fMu_
 125     (
 126         sqr(scalar(1) - exp(-0.0165*(sqrt(k_)*y_/nu())))
 127        *(scalar(1) + 20.5/(Rt_ + SMALL))
 128     ),
 129
 130     nut_
 131     (
 132         IOobject
 133         (
 134             "nut",
 135             runTime_.timeName(),
 136             mesh_,
 137             IOobject::NO_READ,
 138             IOobject::AUTO_WRITE
 139         ),
 140         autoCreateNut("nut", mesh_)
 141     )
 142 {
 143     nut_ = Cmu_*fMu_*sqr(k_)/(epsilon_ + epsilonSmall_);
 144     nut_.correctBoundaryConditions();
 145
 146     printCoeffs();
 147 }
 148
 149
 150 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 151
 152 tmp<volSymmTensorField> LamBremhorstKE::R() const
 153 {
 154     return tmp<volSymmTensorField>
 155     (
 156         new volSymmTensorField
 157         (
 158             IOobject
 159             (
 160                 "R",
 161                 runTime_.timeName(),
 162                 mesh_,
 163                 IOobject::NO_READ,
 164                 IOobject::NO_WRITE
 165             ),
 166             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
 167             k_.boundaryField().types()
 168         )
 169     );
 170 }
 171
 172
 173 tmp<volSymmTensorField> LamBremhorstKE::devReff() const
 174 {
 175     return tmp<volSymmTensorField>
 176     (
 177         new volSymmTensorField
 178         (
 179             IOobject
 180             (
 181                 "devRhoReff",
 182                 runTime_.timeName(),
 183                 mesh_,
 184                 IOobject::NO_READ,
 185                 IOobject::NO_WRITE
 186             ),
 187            -nuEff()*dev(twoSymm(fvc::grad(U_)))
 188         )
 189     );
 190 }
 191
 192
 193 tmp<fvVectorMatrix> LamBremhorstKE::divDevReff(volVectorField& U) const
 194 {
 195     return
 196     (
 197       - fvm::laplacian(nuEff(), U)
 198       - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
 199     );
 200 }
 201
 202
 203 bool LamBremhorstKE::read()
 204 {
 205     if (RASModel::read())
 206     {
 207         Cmu_.readIfPresent(coeffDict());
 208         C1_.readIfPresent(coeffDict());
 209         C2_.readIfPresent(coeffDict());
 210         sigmaEps_.readIfPresent(coeffDict());
 211
 212         return true;
 213     }
 214     else
 215     {
 216         return false;
 217     }
 218 }
 219
 220
 221 void LamBremhorstKE::correct()
 222 {
 223     RASModel::correct();
 224
 225     if (!turbulence_)
 226     {
 227         return;
 228     }
 229
 230     if (mesh_.changing())
 231     {
 232         y_.correct();
 233     }
 234
 235     volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_))));
 236
 237     // Calculate parameters and coefficients for low-Reynolds number model
 238
 239     Rt_ = sqr(k_)/(nu()*epsilon_);
 240     volScalarField Ry = sqrt(k_)*y_/nu();
 241
 242     fMu_ = sqr(scalar(1) - exp(-0.0165*Ry))*(scalar(1) + 20.5/(Rt_ + SMALL));
 243
 244     volScalarField f1 = scalar(1) + pow(0.05/(fMu_ + SMALL), 3);
 245     volScalarField f2 = scalar(1) - exp(-sqr(Rt_));
 246
 247     // Update espsilon and G at the wall
 248     epsilon_.boundaryField().updateCoeffs();
 249
 250     // Dissipation equation
 251
 252     tmp<fvScalarMatrix> epsEqn
 253     (
 254         fvm::ddt(epsilon_)
 255       + fvm::div(phi_, epsilon_)
 256       - fvm::laplacian(DepsilonEff(), epsilon_)
 257      ==
 258         C1_*f1*G*epsilon_/k_
 259       - fvm::Sp(C2_*f2*epsilon_/k_, epsilon_)
 260     );
 261
 262     epsEqn().relax();
 263
 264     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 265
 266     solve(epsEqn);
 267     bound(epsilon_, epsilon0_);
 268
 269
 270     // Turbulent kinetic energy equation
 271
 272     tmp<fvScalarMatrix> kEqn
 273     (
 274         fvm::ddt(k_)
 275       + fvm::div(phi_, k_)
 276       - fvm::laplacian(DkEff(), k_)
 277      ==
 278         G - fvm::Sp(epsilon_/k_, k_)
 279     );
 280
 281     kEqn().relax();
 282     solve(kEqn);
 283     bound(k_, k0_);
 284
 285
 286     // Re-calculate viscosity
 287     nut_ = Cmu_*fMu_*sqr(k_)/epsilon_;
 288     nut_.correctBoundaryConditions();
 289 }
 290
 291
 292 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 293
 294 } // End namespace RASModels
 295 } // End namespace incompressible
 296 } // End namespace Foam
 297
 298 // ************************************************************************* //