src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-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 "radiativeIntensityRay.H"
  28 #include "fvm.H"
  29 #include "fvDOM.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 Foam::label Foam::radiation::radiativeIntensityRay::rayId(0);
  34
  35 const Foam::word
  36 Foam::radiation::radiativeIntensityRay::intensityPrefix("ILambda");
  37
  38
  39 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  40
  41 Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
  42 (
  43     const fvDOM& dom,
  44     const fvMesh& mesh,
  45     const scalar phi,
  46     const scalar theta,
  47     const scalar deltaPhi,
  48     const scalar deltaTheta,
  49     const label nLambda,
  50     const absorptionEmissionModel& absorptionEmission,
  51     const blackBodyEmission& blackBody
  52 )
  53 :
  54     dom_(dom),
  55     mesh_(mesh),
  56     absorptionEmission_(absorptionEmission),
  57     blackBody_(blackBody),
  58     I_
  59     (
  60         IOobject
  61         (
  62             "I" + name(rayId),
  63             mesh_.time().timeName(),
  64             mesh_,
  65             IOobject::NO_READ,
  66             IOobject::NO_WRITE
  67         ),
  68         mesh_,
  69         dimensionedScalar("I", dimMass/pow3(dimTime), 0.0)
  70     ),
  71     Qr_
  72     (
  73         IOobject
  74         (
  75             "Qr" + name(rayId),
  76             mesh_.time().timeName(),
  77             mesh_,
  78             IOobject::NO_READ,
  79             IOobject::NO_WRITE
  80         ),
  81         mesh_,
  82         dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
  83     ),
  84     d_(vector::zero),
  85     dAve_(vector::zero),
  86     theta_(theta),
  87     phi_(phi),
  88     omega_(0.0),
  89     nLambda_(nLambda),
  90     ILambda_(nLambda)
  91 {
  92     scalar sinTheta = Foam::sin(theta);
  93     scalar cosTheta = Foam::cos(theta);
  94     scalar sinPhi = Foam::sin(phi);
  95     scalar cosPhi = Foam::cos(phi);
  96
  97     omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
  98     d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
  99     dAve_ = vector
 100     (
 101         sinPhi
 102        *Foam::sin(0.5*deltaPhi)
 103        *(deltaTheta - Foam::cos(2.0*theta)
 104        *Foam::sin(deltaTheta)),
 105         cosPhi
 106        *Foam::sin(0.5*deltaPhi)
 107        *(deltaTheta - Foam::cos(2.0*theta)
 108        *Foam::sin(deltaTheta)),
 109         0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
 110     );
 111
 112
 113     autoPtr<volScalarField> IDefaultPtr;
 114
 115     forAll(ILambda_, lambdaI)
 116     {
 117         IOobject IHeader
 118         (
 119             intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
 120             mesh_.time().timeName(),
 121             mesh_,
 122             IOobject::MUST_READ,
 123             IOobject::AUTO_WRITE
 124         );
 125
 126         // check if field exists and can be read
 127         if (IHeader.headerOk())
 128         {
 129             ILambda_.set
 130             (
 131                 lambdaI,
 132                 new volScalarField(IHeader, mesh_)
 133             );
 134         }
 135         else
 136         {
 137             // Demand driven load the IDefault field
 138             if (!IDefaultPtr.valid())
 139             {
 140                 IDefaultPtr.reset
 141                 (
 142                     new volScalarField
 143                     (
 144                         IOobject
 145                         (
 146                             "IDefault",
 147                             mesh_.time().timeName(),
 148                             mesh_,
 149                             IOobject::MUST_READ,
 150                             IOobject::NO_WRITE
 151                         ),
 152                         mesh_
 153                     )
 154                 );
 155             }
 156
 157             // Reset the MUST_READ flag
 158             IOobject noReadHeader(IHeader);
 159             noReadHeader.readOpt() = IOobject::NO_READ;
 160
 161             ILambda_.set
 162             (
 163                 lambdaI,
 164                 new volScalarField(noReadHeader, IDefaultPtr())
 165             );
 166         }
 167     }
 168     rayId++;
 169 }
 170
 171
 172 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 173
 174 Foam::radiation::radiativeIntensityRay::~radiativeIntensityRay()
 175 {}
 176
 177
 178 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 179
 180 Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
 181 {
 182     // reset boundary heat flux to zero
 183     //Qr_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
 184     Qr_.boundaryField() = 0.0;
 185
 186     scalar maxResidual = -GREAT;
 187
 188     forAll(ILambda_, lambdaI)
 189     {
 190         const volScalarField& k = dom_.aLambda(lambdaI);
 191
 192         surfaceScalarField Ji = dAve_ & mesh_.Sf();
 193
 194         fvScalarMatrix IiEq
 195         (
 196             fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
 197           + fvm::Sp(k*omega_, ILambda_[lambdaI])
 198          ==
 199             1.0/Foam::mathematicalConstant::pi*omega_
 200            *(
 201                 k*blackBody_.bLambda(lambdaI)
 202               + absorptionEmission_.ECont(lambdaI)
 203             )
 204         );
 205
 206         IiEq.relax();
 207
 208         scalar eqnResidual = solve
 209         (
 210             IiEq,
 211             mesh_.solver("Ii")
 212         ).initialResidual();
 213
 214         maxResidual = max(eqnResidual, maxResidual);
 215     }
 216
 217     return maxResidual;
 218 }
 219
 220
 221 void Foam::radiation::radiativeIntensityRay::addIntensity()
 222 {
 223     I_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
 224
 225     forAll(ILambda_, lambdaI)
 226     {
 227         I_ += absorptionEmission_.addIntensity(lambdaI, ILambda_[lambdaI]);
 228     }
 229 }
 230
 231
 232 // ************************************************************************* //