FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

fvDOM.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "fvDOM.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 
00029 #include <OpenFOAM/mathematicalConstants.H>
00030 #include <radiation/radiationConstants.H>
00031 
00032 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00033 
00034 namespace Foam
00035 {
00036     namespace radiation
00037     {
00038         defineTypeNameAndDebug(fvDOM, 0);
00039 
00040         addToRunTimeSelectionTable
00041         (
00042             radiationModel,
00043             fvDOM,
00044             dictionary
00045         );
00046     }
00047 }
00048 
00049 
00050 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00051 
00052 Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
00053 :
00054     radiationModel(typeName, T),
00055     G_
00056     (
00057         IOobject
00058         (
00059             "G",
00060             mesh_.time().timeName(),
00061             mesh_,
00062             IOobject::NO_READ,
00063             IOobject::AUTO_WRITE
00064         ),
00065         mesh_,
00066         dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
00067     ),
00068     Qr_
00069     (
00070         IOobject
00071         (
00072             "Qr",
00073             mesh_.time().timeName(),
00074             mesh_,
00075             IOobject::NO_READ,
00076             IOobject::AUTO_WRITE
00077         ),
00078         mesh_,
00079         dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
00080     ),
00081     Qem_
00082     (
00083         IOobject
00084         (
00085             "Qem",
00086             mesh_.time().timeName(),
00087             mesh_,
00088             IOobject::NO_READ,
00089             IOobject::NO_WRITE
00090         ),
00091         mesh_,
00092         dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
00093     ),
00094     Qin_
00095     (
00096         IOobject
00097         (
00098             "Qin",
00099             mesh_.time().timeName(),
00100             mesh_,
00101             IOobject::NO_READ,
00102             IOobject::NO_WRITE
00103         ),
00104         mesh_,
00105         dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
00106     ),
00107     a_
00108     (
00109         IOobject
00110         (
00111             "a",
00112             mesh_.time().timeName(),
00113             mesh_,
00114             IOobject::NO_READ,
00115             IOobject::NO_WRITE
00116         ),
00117         mesh_,
00118         dimensionedScalar("a", dimless/dimLength, 0.0)
00119     ),
00120     e_
00121     (
00122         IOobject
00123         (
00124             "e",
00125             mesh_.time().timeName(),
00126             mesh_,
00127             IOobject::NO_READ,
00128             IOobject::NO_WRITE
00129         ),
00130         mesh_,
00131         dimensionedScalar("a", dimless/dimLength, 0.0)
00132     ),
00133     E_
00134     (
00135         IOobject
00136         (
00137             "E",
00138             mesh_.time().timeName(),
00139             mesh_,
00140             IOobject::NO_READ,
00141             IOobject::NO_WRITE
00142         ),
00143         mesh_,
00144         dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
00145     ),
00146     nTheta_(readLabel(coeffs_.lookup("nTheta"))),
00147     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
00148     nRay_(0),
00149     nLambda_(absorptionEmission_->nBands()),
00150     aLambda_(nLambda_),
00151     blackBody_(nLambda_, T),
00152     IRay_(0),
00153     convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
00154     maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50))
00155 {
00156     if (mesh_.nSolutionD() == 3)    //3D
00157     {
00158         nRay_ = 4*nPhi_*nTheta_;
00159         IRay_.setSize(nRay_);
00160         scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
00161         scalar deltaTheta = mathematicalConstant::pi/nTheta_;
00162         label i = 0;
00163         for (label n = 1; n <= nTheta_; n++)
00164         {
00165             for (label m = 1; m <= 4*nPhi_; m++)
00166             {
00167                 scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
00168                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
00169                 IRay_.set
00170                 (
00171                     i,
00172                     new radiativeIntensityRay
00173                     (
00174                         *this,
00175                         mesh_,
00176                         phii,
00177                         thetai,
00178                         deltaPhi,
00179                         deltaTheta,
00180                         nLambda_,
00181                         absorptionEmission_,
00182                         blackBody_,
00183                         i
00184                     )
00185                 );
00186                 i++;
00187             }
00188         }
00189     }
00190     else
00191     {
00192         if (mesh_.nSolutionD() == 2)    //2D (X & Y)
00193         {
00194             scalar thetai = mathematicalConstant::piByTwo;
00195             scalar deltaTheta = mathematicalConstant::pi;
00196             nRay_ = 4*nPhi_;
00197             IRay_.setSize(nRay_);
00198             scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
00199             label i = 0;
00200             for (label m = 1; m <= 4*nPhi_; m++)
00201             {
00202                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
00203                 IRay_.set
00204                 (
00205                     i,
00206                     new radiativeIntensityRay
00207                     (
00208                         *this,
00209                         mesh_,
00210                         phii,
00211                         thetai,
00212                         deltaPhi,
00213                         deltaTheta,
00214                         nLambda_,
00215                         absorptionEmission_,
00216                         blackBody_,
00217                         i
00218                     )
00219                 );
00220                 i++;
00221             }
00222         }
00223         else    //1D (X)
00224         {
00225             scalar thetai = mathematicalConstant::piByTwo;
00226             scalar deltaTheta = mathematicalConstant::pi;
00227             nRay_ = 2;
00228             IRay_.setSize(nRay_);
00229             scalar deltaPhi = mathematicalConstant::pi;
00230             label i = 0;
00231             for (label m = 1; m <= 2; m++)
00232             {
00233                 scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
00234                 IRay_.set
00235                 (
00236                     i,
00237                     new radiativeIntensityRay
00238                     (
00239                         *this,
00240                         mesh_,
00241                         phii,
00242                         thetai,
00243                         deltaPhi,
00244                         deltaTheta,
00245                         nLambda_,
00246                         absorptionEmission_,
00247                         blackBody_,
00248                         i
00249                     )
00250                 );
00251                 i++;
00252             }
00253 
00254         }
00255     }
00256 
00257 
00258     // Construct absorption field for each wavelength
00259     forAll(aLambda_, lambdaI)
00260     {
00261         aLambda_.set
00262         (
00263             lambdaI,
00264             new volScalarField
00265             (
00266                 IOobject
00267                 (
00268                     "aLambda_" + Foam::name(lambdaI) ,
00269                     mesh_.time().timeName(),
00270                     mesh_,
00271                     IOobject::NO_READ,
00272                     IOobject::NO_WRITE
00273                 ),
00274                 a_
00275             )
00276         );
00277     }
00278 
00279     Info<< "fvDOM : Allocated " << IRay_.size()
00280         << " rays with average orientation:" << nl;
00281     forAll (IRay_, i)
00282     {
00283         Info<< '\t' << IRay_[i].I().name()
00284             << '\t' << IRay_[i].dAve() << nl;
00285     }
00286     Info<< endl;
00287 }
00288 
00289 
00290 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00291 
00292 Foam::radiation::fvDOM::~fvDOM()
00293 {}
00294 
00295 
00296 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00297 
00298 bool Foam::radiation::fvDOM::read()
00299 {
00300     if (radiationModel::read())
00301     {
00302 //      Only reading solution parameters - not changing ray geometry
00303 
00304         coeffs_.readIfPresent("convergence", convergence_);
00305         coeffs_.readIfPresent("maxIter", maxIter_);
00306 
00307         return true;
00308     }
00309     else
00310     {
00311         return false;
00312     }
00313 }
00314 
00315 
00316 void Foam::radiation::fvDOM::calculate()
00317 {
00318     absorptionEmission_->correct(a_, aLambda_);
00319 
00320     updateBlackBodyEmission();
00321 
00322     scalar maxResidual = 0.0;
00323     label radIter = 0;
00324     do
00325     {
00326         radIter++;
00327         forAll(IRay_, rayI)
00328         {
00329             maxResidual = 0.0;
00330             scalar maxBandResidual = IRay_[rayI].correct();
00331             maxResidual = max(maxBandResidual, maxResidual);
00332         }
00333 
00334         Info << "Radiation solver iter: " << radIter << endl;
00335 
00336     } while(maxResidual > convergence_ && radIter < maxIter_);
00337 
00338     updateG();
00339 }
00340 
00341 
00342 Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
00343 {
00344     return tmp<volScalarField>
00345     (
00346         new volScalarField
00347         (
00348             IOobject
00349             (
00350                 "Rp",
00351                 mesh_.time().timeName(),
00352                 mesh_,
00353                 IOobject::NO_READ,
00354                 IOobject::NO_WRITE,
00355                 false
00356             ),
00357             4.0*a_*radiation::sigmaSB //absorptionEmission_->a()
00358         )
00359     );
00360 }
00361 
00362 
00363 Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
00364 Foam::radiation::fvDOM::Ru() const
00365 {
00366 
00367     const DimensionedField<scalar, volMesh>& G =
00368         G_.dimensionedInternalField();
00369     const DimensionedField<scalar, volMesh> E =
00370         absorptionEmission_->ECont()().dimensionedInternalField();
00371     const DimensionedField<scalar, volMesh> a =
00372         a_.dimensionedInternalField(); //absorptionEmission_->aCont()()
00373 
00374     return  a*G - E;
00375 }
00376 
00377 
00378 void Foam::radiation::fvDOM::updateBlackBodyEmission()
00379 {
00380     for (label j=0; j < nLambda_; j++)
00381     {
00382         blackBody_.correct(j, absorptionEmission_->bands(j));
00383     }
00384 }
00385 
00386 
00387 void Foam::radiation::fvDOM::updateG()
00388 {
00389     G_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00390     Qr_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00391     Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00392     Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
00393 
00394     forAll(IRay_, rayI)
00395     {
00396         IRay_[rayI].addIntensity();
00397         G_ += IRay_[rayI].I()*IRay_[rayI].omega();
00398         Qr_.boundaryField() += IRay_[rayI].Qr().boundaryField();
00399         Qem_.boundaryField() += IRay_[rayI].Qem().boundaryField();
00400         Qin_.boundaryField() += IRay_[rayI].Qin().boundaryField();
00401     }
00402 }
00403 
00404 
00405 void Foam::radiation::fvDOM::setRayIdLambdaId
00406 (
00407     const word& name,
00408     label& rayId,
00409     label& lambdaId
00410 ) const
00411 {
00412     // assuming name is in the form: CHARS_rayId_lambdaId
00413     size_type i1 = name.find_first_of("_");
00414     size_type i2 = name.find_last_of("_");
00415 
00416     rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
00417     lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
00418 }
00419 
00420 
00421 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines