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