applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.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 "SCOPEXiEq.H"
  28 #include "addToRunTimeSelectionTable.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34 namespace XiEqModels
  35 {
  36     defineTypeNameAndDebug(SCOPEXiEq, 0);
  37     addToRunTimeSelectionTable(XiEqModel, SCOPEXiEq, dictionary);
  38 };
  39 };
  40
  41
  42 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  43
  44 Foam::XiEqModels::SCOPEXiEq::SCOPEXiEq
  45 (
  46     const dictionary& XiEqProperties,
  47     const hhuCombustionThermo& thermo,
  48     const compressible::RASModel& turbulence,
  49     const volScalarField& Su
  50 )
  51 :
  52     XiEqModel(XiEqProperties, thermo, turbulence, Su),
  53     XiEqCoef(readScalar(XiEqModelCoeffs_.lookup("XiEqCoef"))),
  54     XiEqExp(readScalar(XiEqModelCoeffs_.lookup("XiEqExp"))),
  55     lCoef(readScalar(XiEqModelCoeffs_.lookup("lCoef"))),
  56     SuMin(0.01*Su.average()),
  57     MaModel
  58     (
  59         IOdictionary
  60         (
  61             IOobject
  62             (
  63                 "combustionProperties",
  64                 Su.mesh().time().constant(),
  65                 Su.mesh(),
  66                 IOobject::MUST_READ
  67             )
  68         ),
  69         thermo
  70     )
  71 {}
  72
  73
  74 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
  75
  76 Foam::XiEqModels::SCOPEXiEq::~SCOPEXiEq()
  77 {}
  78
  79
  80 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
  81
  82 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
  83 {
  84     const volScalarField& k = turbulence_.k();
  85     const volScalarField& epsilon = turbulence_.epsilon();
  86
  87     volScalarField up = sqrt((2.0/3.0)*k);
  88     volScalarField l = (lCoef*sqrt(3.0/2.0))*up*k/epsilon;
  89     volScalarField Rl = up*l*thermo_.rhou()/thermo_.muu();
  90
  91     volScalarField upBySu = up/(Su_ + SuMin);
  92     volScalarField K = 0.157*upBySu/sqrt(Rl);
  93     volScalarField Ma = MaModel.Ma();
  94
  95     tmp<volScalarField> tXiEq
  96     (
  97         new volScalarField
  98         (
  99             IOobject
 100             (
 101                 "XiEq",
 102                 epsilon.time().timeName(),
 103                 epsilon.db(),
 104                 IOobject::NO_READ,
 105                 IOobject::NO_WRITE
 106             ),
 107             epsilon.mesh(),
 108             dimensionedScalar("XiEq", dimless, 0.0)
 109         )
 110     );
 111     volScalarField& xieq = tXiEq();
 112
 113     forAll(xieq, celli)
 114     {
 115         if (Ma[celli] > 0.01)
 116         {
 117             xieq[celli] =
 118                 XiEqCoef*pow(K[celli]*Ma[celli], -XiEqExp)*upBySu[celli];
 119         }
 120     }
 121
 122     forAll(xieq.boundaryField(), patchi)
 123     {
 124         scalarField& xieqp = xieq.boundaryField()[patchi];
 125         const scalarField& Kp = K.boundaryField()[patchi];
 126         const scalarField& Map = Ma.boundaryField()[patchi];
 127         const scalarField& upBySup = upBySu.boundaryField()[patchi];
 128
 129         forAll(xieqp, facei)
 130         {
 131             if (Ma[facei] > 0.01)
 132             {
 133                 xieqp[facei] =
 134                     XiEqCoef*pow(Kp[facei]*Map[facei], -XiEqExp)*upBySup[facei];
 135             }
 136         }
 137     }
 138
 139     return tXiEq;
 140 }
 141
 142
 143 bool Foam::XiEqModels::SCOPEXiEq::read(const dictionary& XiEqProperties)
 144 {
 145     XiEqModel::read(XiEqProperties);
 146
 147     XiEqModelCoeffs_.lookup("XiEqCoef") >> XiEqCoef;
 148     XiEqModelCoeffs_.lookup("XiEqExp") >> XiEqExp;
 149     XiEqModelCoeffs_.lookup("lCoef") >> lCoef;
 150
 151     return true;
 152 }
 153
 154
 155 // ************************************************************************* //