src/turbulenceModels/RAS/incompressible/RASModel/RASModel.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 "RASModel.H"
  28 #include "wallFvPatch.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34 namespace incompressible
  35 {
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 defineTypeNameAndDebug(RASModel, 0);
  40 defineRunTimeSelectionTable(RASModel, dictionary);
  41
  42 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
  43
  44 void RASModel::printCoeffs()
  45 {
  46     if (printCoeffs_)
  47     {
  48         Info<< type() << "Coeffs" << coeffDict_ << endl;;
  49     }
  50 }
  51
  52
  53 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  54
  55 RASModel::RASModel
  56 (
  57     const word& type,
  58     const volVectorField& U,
  59     const surfaceScalarField& phi,
  60     transportModel& lamTransportModel
  61 )
  62 :
  63     IOdictionary
  64     (
  65         IOobject
  66         (
  67             "RASProperties",
  68             U.time().constant(),
  69             U.db(),
  70             IOobject::MUST_READ,
  71             IOobject::NO_WRITE
  72         )
  73     ),
  74
  75     runTime_(U.time()),
  76     mesh_(U.mesh()),
  77
  78     U_(U),
  79     phi_(phi),
  80     transportModel_(lamTransportModel),
  81
  82     turbulence_(lookup("turbulence")),
  83     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
  84     coeffDict_(subDict(type + "Coeffs")),
  85
  86     kappa_
  87     (
  88         dimensioned<scalar>::lookupOrAddToDict
  89         (
  90             "kappa",
  91             subDict("wallFunctionCoeffs"),
  92             0.4187
  93         )
  94     ),
  95     E_
  96     (
  97         dimensioned<scalar>::lookupOrAddToDict
  98         (
  99             "E",
 100             subDict("wallFunctionCoeffs"),
 101             9.0
 102         )
 103     ),
 104
 105     yPlusLam_(yPlusLam(kappa_.value(), E_.value())),
 106
 107     k0_("k0", dimVelocity*dimVelocity, SMALL),
 108     epsilon0_("epsilon", k0_.dimensions()/dimTime, SMALL),
 109     epsilonSmall_("epsilonSmall", epsilon0_.dimensions(), SMALL),
 110
 111     y_(mesh_)
 112 {}
 113
 114
 115 RASModel::~RASModel()
 116 {}
 117
 118
 119 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 120
 121 scalar RASModel::yPlusLam(const scalar kappa, const scalar E)
 122 {
 123     scalar ypl = 11.0;
 124
 125     for (int i=0; i<10; i++)
 126     {
 127         ypl = log(E*ypl)/kappa;
 128     }
 129
 130     return ypl;
 131 }
 132
 133
 134 tmp<scalarField> RASModel::yPlus(const label patchNo) const
 135 {
 136     const fvPatch& curPatch = mesh_.boundary()[patchNo];
 137
 138     tmp<scalarField> tYp(new scalarField(curPatch.size()));
 139     scalarField& Yp = tYp();
 140
 141     if (typeid(curPatch) == typeid(wallFvPatch))
 142     {
 143         scalar Cmu(readScalar(coeffDict_.lookup("Cmu")));
 144
 145         Yp = pow(Cmu, 0.25)*y_[patchNo]
 146             *sqrt(k()().boundaryField()[patchNo].patchInternalField())
 147             /nu().boundaryField()[patchNo];
 148     }
 149     else
 150     {
 151         WarningIn
 152         (
 153             "tmp<scalarField> RASModel::yPlus(const label patchNo)"
 154         )   << "const : " << endl
 155             << "Patch " << patchNo << " is not a wall.  Returning blank field"
 156             << endl;
 157
 158         Yp.setSize(0);
 159     }
 160
 161     return tYp;
 162 }
 163
 164
 165 void RASModel::correct()
 166 {
 167     if (mesh_.changing())
 168     {
 169         y_.correct();
 170     }
 171 }
 172
 173
 174 bool RASModel::read()
 175 {
 176     if (regIOobject::read())
 177     {
 178         lookup("turbulence") >> turbulence_;
 179         coeffDict_ = subDict(type() + "Coeffs");
 180
 181         kappa_.readIfPresent(subDict("wallFunctionCoeffs"));
 182         E_.readIfPresent(subDict("wallFunctionCoeffs"));
 183
 184         yPlusLam_ = yPlusLam(kappa_.value(), E_.value());
 185
 186         k0_.readIfPresent(*this);
 187         epsilon0_.readIfPresent(*this);
 188         epsilonSmall_.readIfPresent(*this);
 189
 190         return true;
 191     }
 192     else
 193     {
 194         return false;
 195     }
 196 }
 197
 198
 199 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 200
 201 } // End namespace incompressible
 202 } // End namespace Foam
 203
 204 // ************************************************************************* //