src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.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 "fvDOM.H"
  28 #include "addToRunTimeSelectionTable.H"
  29
  30 #include "mathematicalConstants.H"
  31 #include "radiationConstants.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37     namespace radiation
  38     {
  39         defineTypeNameAndDebug(fvDOM, 0);
  40
  41         addToRunTimeSelectionTable
  42         (
  43             radiationModel,
  44             fvDOM,
  45             dictionary
  46         );
  47     }
  48 }
  49
  50
  51 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  52
  53 Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
  54 :
  55     radiationModel(typeName, T),
  56     G_
  57     (
  58         IOobject
  59         (
  60             "G",
  61             mesh_.time().timeName(),
  62             mesh_,
  63             IOobject::NO_READ,
  64             IOobject::AUTO_WRITE
  65         ),
  66         mesh_,
  67         dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
  68     ),
  69     Qr_
  70     (
  71         IOobject
  72         (
  73             "Qr",
  74             mesh_.time().timeName(),
  75             mesh_,
  76             IOobject::NO_READ,
  77             IOobject::AUTO_WRITE
  78         ),
  79         mesh_,
  80         dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
  81     ),
  82     a_
  83     (
  84         IOobject
  85         (
  86             "a",
  87             mesh_.time().timeName(),
  88             mesh_,
  89             IOobject::NO_READ,
  90             IOobject::AUTO_WRITE
  91         ),
  92         mesh_,
  93         dimensionedScalar("a", dimless/dimLength, 0.0)
  94     ),
  95     e_
  96     (
  97         IOobject
  98         (
  99             "e",
 100             mesh_.time().timeName(),
 101             mesh_,
 102             IOobject::NO_READ,
 103             IOobject::NO_WRITE
 104         ),
 105         mesh_,
 106         dimensionedScalar("a", dimless/dimLength, 0.0)
 107     ),
 108     E_
 109     (
 110         IOobject
 111         (
 112             "E",
 113             mesh_.time().timeName(),
 114             mesh_,
 115             IOobject::NO_READ,
 116             IOobject::NO_WRITE
 117         ),
 118         mesh_,
 119         dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
 120     ),
 121     nTheta_(readLabel(coeffs_.lookup("nTheta"))),
 122     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
 123     nRay_(0),
 124     nLambda_(absorptionEmission_->nBands()),
 125     aLambda_(nLambda_),
 126     blackBody_(nLambda_, T),
 127     IRay_(0),
 128     convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
 129     maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50))
 130 {
 131     if (mesh_.nSolutionD() == 3)    //3D
 132     {
 133         nRay_ = 4*nPhi_*nTheta_;
 134         IRay_.setSize(nRay_);
 135         scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
 136         scalar deltaTheta = mathematicalConstant::pi/nTheta_;
 137         label i = 0;
 138         for (label n = 1; n <= nTheta_; n++)
 139         {
 140             for (label m = 1; m <= 4*nPhi_; m++)
 141             {
 142                 scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
 143                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
 144                 IRay_.set
 145                 (
 146                     i,
 147                     new radiativeIntensityRay
 148                     (
 149                         *this,
 150                         mesh_,
 151                         phii,
 152                         thetai,
 153                         deltaPhi,
 154                         deltaTheta,
 155                         nLambda_,
 156                         absorptionEmission_,
 157                         blackBody_
 158                     )
 159                 );
 160                 i++;
 161             }
 162         }
 163     }
 164     else
 165     {
 166         if (mesh_.nSolutionD() == 2)    //2D (X & Y)
 167         {
 168             scalar thetai = mathematicalConstant::piByTwo;
 169             scalar deltaTheta = mathematicalConstant::pi;
 170             nRay_ = 4*nPhi_;
 171             IRay_.setSize(nRay_);
 172             scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
 173             label i = 0;
 174             for (label m = 1; m <= 4*nPhi_; m++)
 175             {
 176                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
 177                 IRay_.set
 178                 (
 179                     i,
 180                     new radiativeIntensityRay
 181                     (
 182                         *this,
 183                         mesh_,
 184                         phii,
 185                         thetai,
 186                         deltaPhi,
 187                         deltaTheta,
 188                         nLambda_,
 189                         absorptionEmission_,
 190                         blackBody_
 191                     )
 192                 );
 193                 i++;
 194             }
 195         }
 196         else    //1D (X)
 197         {
 198             scalar thetai = mathematicalConstant::piByTwo;
 199             scalar deltaTheta = mathematicalConstant::pi;
 200             nRay_ = 2;
 201             IRay_.setSize(nRay_);
 202             scalar deltaPhi = mathematicalConstant::pi;
 203             label i = 0;
 204             for (label m = 1; m <= 2; m++)
 205             {
 206                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
 207                 IRay_.set
 208                 (
 209                     i,
 210                     new radiativeIntensityRay
 211                     (
 212                         *this,
 213                         mesh_,
 214                         phii,
 215                         thetai,
 216                         deltaPhi,
 217                         deltaTheta,
 218                         nLambda_,
 219                         absorptionEmission_,
 220                         blackBody_
 221                     )
 222                 );
 223                 i++;
 224             }
 225
 226         }
 227     }
 228
 229
 230     // Construct absorption field for each wavelength
 231     forAll(aLambda_, lambdaI)
 232     {
 233         aLambda_.set
 234         (
 235             lambdaI,
 236             new volScalarField
 237             (
 238                 IOobject
 239                 (
 240                     "aLambda_" + Foam::name(lambdaI) ,
 241                     mesh_.time().timeName(),
 242                     mesh_,
 243                     IOobject::NO_READ,
 244                     IOobject::NO_WRITE
 245                 ),
 246                 a_
 247             )
 248         );
 249     }
 250
 251     Info<< "fvDOM : Allocated " << IRay_.size()
 252         << " rays with average orientation:" << nl;
 253     forAll (IRay_, i)
 254     {
 255         Info<< '\t' << IRay_[i].I().name()
 256             << '\t' << IRay_[i].dAve() << nl;
 257     }
 258     Info<< endl;
 259 }
 260
 261
 262 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 263
 264 Foam::radiation::fvDOM::~fvDOM()
 265 {}
 266
 267
 268 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 269
 270 bool Foam::radiation::fvDOM::read()
 271 {
 272     if (radiationModel::read())
 273     {
 274 //      Only reading solution parameters - not changing ray geometry
 275
 276         coeffs_.readIfPresent("convergence", convergence_);
 277         coeffs_.readIfPresent("maxIter", maxIter_);
 278
 279         return true;
 280     }
 281     else
 282     {
 283         return false;
 284     }
 285 }
 286
 287
 288 void Foam::radiation::fvDOM::calculate()
 289 {
 290     absorptionEmission_->correct(a_, aLambda_);
 291
 292     updateBlackBodyEmission();
 293
 294     scalar maxResidual = 0.0;
 295     label radIter = 0;
 296     do
 297     {
 298         radIter++;
 299         forAll(IRay_, rayI)
 300         {
 301             maxResidual = 0.0;
 302             scalar maxBandResidual = IRay_[rayI].correct();
 303             maxResidual = max(maxBandResidual, maxResidual);
 304         }
 305
 306         Info << "Radiation solver iter: " << radIter << endl;
 307
 308     } while(maxResidual > convergence_ && radIter < maxIter_);
 309
 310     updateG();
 311 }
 312
 313
 314 Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
 315 {
 316     return tmp<volScalarField>
 317     (
 318         new volScalarField
 319         (
 320             IOobject
 321             (
 322                 "Rp",
 323                 mesh_.time().timeName(),
 324                 mesh_,
 325                 IOobject::NO_READ,
 326                 IOobject::NO_WRITE,
 327                 false
 328             ),
 329             4.0*a_*radiation::sigmaSB //absorptionEmission_->a()
 330         )
 331     );
 332 }
 333
 334
 335 Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
 336 Foam::radiation::fvDOM::Ru() const
 337 {
 338
 339     const DimensionedField<scalar, volMesh>& G =
 340         G_.dimensionedInternalField();
 341     const DimensionedField<scalar, volMesh> E =
 342         absorptionEmission_->ECont()().dimensionedInternalField();
 343     const DimensionedField<scalar, volMesh> a =
 344         a_.dimensionedInternalField(); //absorptionEmission_->aCont()()
 345
 346     return  a*G - E;
 347 }
 348
 349
 350 void Foam::radiation::fvDOM::updateBlackBodyEmission()
 351 {
 352     for (label j=0; j < nLambda_; j++)
 353     {
 354         blackBody_.correct(j, absorptionEmission_->bands(j));
 355     }
 356 }
 357
 358
 359 void Foam::radiation::fvDOM::updateG()
 360 {
 361     G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
 362     Qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
 363
 364     forAll(IRay_, rayI)
 365     {
 366         IRay_[rayI].addIntensity();
 367         G_ += IRay_[rayI].I()*IRay_[rayI].omega();
 368         //Qr_ += IRay_[rayI].Qr();
 369         Qr_.boundaryField() += IRay_[rayI].Qr().boundaryField();
 370     }
 371 }
 372
 373
 374 void Foam::radiation::fvDOM::setRayIdLambdaId
 375 (
 376     const word& name,
 377     label& rayId,
 378     label& lambdaId
 379 ) const
 380 {
 381     // assuming name is in the form: CHARS_rayId_lambdaId
 382     size_type i1 = name.find_first_of("_");
 383     size_type i2 = name.find_last_of("_");
 384
 385     rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
 386     lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
 387 }
 388
 389
 390 // ************************************************************************* //