ENH: Adding collimate beam source to fvDOM and multiBandZoneAbsorption

src/thermophysicalModels/radiation/Make/files

@@ -41,6 +41,7 @@ submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionE
 submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
 submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
 submodels/absorptionEmissionModel/multiBandAbsorptionEmission/multiBandAbsorptionEmission.C
+submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.C
 submodels/boundaryRadiationProperties/boundaryRadiationProperties.C
 submodels/boundaryRadiationProperties/boundaryRadiationPropertiesPatch.C
 

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -33,6 +33,7 @@ License
 
 #include "fvDOM.H"
 #include "constants.H"
+#include "unitConversion.H"
 
 using namespace Foam::constant;
 using namespace Foam::constant::mathematical;
@@ -98,6 +99,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 
         fvPatchScalarField::operator=(refValue());
     }
+
 }
 
 
@@ -178,6 +180,13 @@ updateCoeffs()
 
     const scalarField& emissivity = temissivity();
 
+    const tmp<scalarField> ttransmissivity
+    (
+        boundaryRadiation.transmissivity(patch().index())
+    );
+
+    const scalarField& transmissivity = ttransmissivity();
+
     scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
     scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
@@ -211,6 +220,39 @@ updateCoeffs()
         }
     }
 
+    scalarField Iexternal(this->size(), 0.0);
+
+    if (dom.useExternalBeam())
+    {
+        const vector sunDir = dom.solarCalc().direction();
+        const scalar directSolarRad = dom.solarCalc().directSolarRad();
+
+        //label nRaysBeam = dom.nRaysBeam();
+        label SunRayId(-1);
+        scalar maxSunRay = -GREAT;
+
+        // Looking for the ray closest to the Sun direction
+        for (label rayI=0; rayI < dom.nRay(); rayI++)
+        {
+            const vector& iD = dom.IRay(rayI).d();
+            scalar dir = sunDir & iD;
+            if (dir > maxSunRay)
+            {
+                maxSunRay = dir;
+                SunRayId = rayI;
+            }
+        }
+
+        if (rayId == SunRayId)
+        {
+            const scalarField nAve(n & dom.IRay(rayId).dAve());
+            forAll(Iexternal, faceI)
+            {
+                Iexternal[faceI] = directSolarRad/mag(dom.IRay(rayId).dAve());
+            }
+        }
+    }
+
     forAll(Iw, faceI)
     {
         if ((-n[faceI] & myRayId) > 0.0)
@@ -219,7 +261,8 @@ updateCoeffs()
             refGrad()[faceI] = 0.0;
             valueFraction()[faceI] = 1.0;
             refValue()[faceI] =
-                (
+                Iexternal[faceI]*transmissivity[faceI]
+              + (
                     Ir[faceI]*(scalar(1) - emissivity[faceI])
                   + emissivity[faceI]*physicoChemical::sigma.value()
                   * pow4(Tp[faceI])

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C

@@ -78,7 +78,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF)
 {
-    if (dict.found("value"))
+    if (dict.found("refValue"))
     {
         fvPatchScalarField::operator=
         (
@@ -176,9 +176,14 @@ updateCoeffs()
     (
         boundaryRadiation.emissivity(patch().index(), lambdaId)
     );
-
     const scalarField& emissivity = temissivity();
 
+    const tmp<scalarField> ttransmissivity
+    (
+        boundaryRadiation.transmissivity(patch().index(), lambdaId)
+    );
+    const scalarField& transmissivity = ttransmissivity();
+
     scalarField& qem = ray.qem().boundaryFieldRef()[patchi];
     scalarField& qin = ray.qin().boundaryFieldRef()[patchi];
 
@@ -232,6 +237,41 @@ updateCoeffs()
         }
     }
 
+    scalarField Iexternal(this->size(), 0.0);
+
+    if (dom.useExternalBeam())
+    {
+        const vector sunDir = dom.solarCalc().direction();
+        const scalar directSolarRad =
+            dom.solarCalc().directSolarRad()
+           *dom.spectralDistribution()[lambdaId];
+
+        //label nRaysBeam = dom.nRaysBeam();
+        label SunRayId(-1);
+        scalar maxSunRay = -GREAT;
+
+        // Looking for the ray closest to the Sun direction
+        for (label rayI=0; rayI < dom.nRay(); rayI++)
+        {
+            const vector& iD = dom.IRay(rayI).d();
+            scalar dir = sunDir & iD;
+            if (dir > maxSunRay)
+            {
+                maxSunRay = dir;
+                SunRayId = rayI;
+            }
+        }
+
+        if (rayId == SunRayId)
+        {
+            const scalarField nAve(n & dom.IRay(rayId).dAve());
+            forAll(Iexternal, faceI)
+            {
+                Iexternal[faceI] = directSolarRad/mag(dom.IRay(rayId).dAve());
+            }
+        }
+    }
+
     forAll(Iw, facei)
     {
         const vector& d = dom.IRay(rayId).d();
@@ -242,7 +282,8 @@ updateCoeffs()
             refGrad()[facei] = 0.0;
             valueFraction()[facei] = 1.0;
             refValue()[facei] =
-                (
+                Iexternal[facei]*transmissivity[facei]
+              + (
                     Ir[facei]*(1.0 - emissivity[facei])
                   + emissivity[facei]*Eb[facei]
                 )/pi;

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C

@@ -29,6 +29,7 @@ License
 #include "absorptionEmissionModel.H"
 #include "scatterModel.H"
 #include "constants.H"
+#include "unitConversion.H"
 #include "fvm.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -49,22 +50,125 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::radiation::fvDOM::rotateInitialRays(const vector& sunDir)
+{
+    // Rotate Y spherical cordinates to Sun direction.
+    // Solid angles on the equator are better fit for planar radiation
+    const tensor coordRot = rotationTensor(vector(0, 1, 0), sunDir);
+
+    forAll(IRay_, rayId)
+    {
+        IRay_[rayId].dAve() = coordRot & IRay_[rayId].dAve();
+        IRay_[rayId].d() = coordRot & IRay_[rayId].d();
+    }
+}
+
+
+void Foam::radiation::fvDOM:: alignClosestRayToSun(const vector& sunDir)
+{
+    label SunRayId(-1);
+    scalar maxSunRay = -GREAT;
+
+    // Looking for the ray closest to the Sun direction
+    forAll(IRay_, rayId)
+    {
+        const vector& iD = IRay_[rayId].d();
+        scalar dir = sunDir & iD;
+        if (dir > maxSunRay)
+        {
+            maxSunRay = dir;
+            SunRayId = rayId;
+        }
+    }
+
+    // Second rotation to align colimated radiation with the closest ray
+    const tensor coordRot = rotationTensor(IRay_[SunRayId].d(), sunDir);
+
+    forAll(IRay_, rayId)
+    {
+        IRay_[rayId].dAve() = coordRot & IRay_[rayId].dAve();
+        IRay_[rayId].d() = coordRot & IRay_[rayId].d();
+    }
+
+    Info << "Sun direction : " << sunDir << nl << endl;
+    Info << "Sun ray ID : " << SunRayId << nl << endl;
+}
+
+
+void Foam::radiation::fvDOM::updateRaysDir()
+{
+    solarCalculator_->correctSunDirection();
+    const vector sunDir = solarCalculator_->direction();
+
+    // First iteration
+    if (updateTimeIndex_ == 0)
+    {
+        rotateInitialRays(sunDir);
+        alignClosestRayToSun(sunDir);
+    }
+    else if (updateTimeIndex_ > 0)
+    {
+        alignClosestRayToSun(sunDir);
+    }
+}
+
+
 void Foam::radiation::fvDOM::initialise()
 {
+    coeffs_.readIfPresent("useExternalBeam", useExternalBeam_);
+
+    if (useExternalBeam_)
+    {
+        coeffs_.readEntry("spectralDistribution", spectralDistribution_);
+
+        spectralDistribution_ =
+            spectralDistribution_/sum(spectralDistribution_);
+
+        const dictionary& solarDict = this->subDict("solarCalculatorCoeffs");
+        solarCalculator_.reset(new solarCalculator(solarDict, mesh_));
+
+        if (mesh_.nSolutionD() != 3)
+        {
+            FatalErrorInFunction
+                << "External beam model only available in 3D meshes "
+                << abort(FatalError);
+        }
+
+        if (solarCalculator_->diffuseSolarRad() > 0)
+        {
+            FatalErrorInFunction
+                << "External beam model does not support Diffuse "
+                << "Solar Radiation. Set diffuseSolarRad to zero"
+                << abort(FatalError);
+        }
+        if (spectralDistribution_.size() != nLambda_)
+        {
+            FatalErrorInFunction
+                << "The epectral energy distribution has different bands "
+                << "than the absoprtivity model "
+                << abort(FatalError);
+        }
+    }
+
     // 3D
     if (mesh_.nSolutionD() == 3)
     {
         nRay_ = 4*nPhi_*nTheta_;
+
         IRay_.setSize(nRay_);
-        const scalar deltaPhi = pi/(2.0*nPhi_);
+
+        const scalar deltaPhi = pi/(2*nPhi_);
         const scalar deltaTheta = pi/nTheta_;
+
         label i = 0;
+
         for (label n = 1; n <= nTheta_; n++)
         {
             for (label m = 1; m <= 4*nPhi_; m++)
             {
-                const scalar thetai = (2*n - 1)*deltaTheta/2.0;
+                scalar thetai = (2*n - 1)*deltaTheta/2.0;
-                const scalar phii = (2*m - 1)*deltaPhi/2.0;
+                scalar phii = (2*m - 1)*deltaPhi/2.0;
+
                 IRay_.set
                 (
                     i,
@@ -176,22 +280,41 @@ void Foam::radiation::fvDOM::initialise()
     Info<< "fvDOM : Allocated " << IRay_.size()
         << " rays with average orientation:" << nl;
 
+    if (useExternalBeam_)
+    {
+        // Rotate rays for Sun direction
+        updateRaysDir();
+    }
+
+    scalar totalOmega = 0;
     forAll(IRay_, rayId)
     {
         if (omegaMax_ <  IRay_[rayId].omega())
         {
             omegaMax_ = IRay_[rayId].omega();
         }
+        totalOmega += IRay_[rayId].omega();
         Info<< '\t' << IRay_[rayId].I().name() << " : " << "dAve : "
-            << '\t' << IRay_[rayId].dAve() << nl;
+            << '\t' << IRay_[rayId].dAve() << " : " << "omega : "
+            << '\t' << IRay_[rayId].omega() << " : " << "d : "
+            << '\t' << IRay_[rayId].d() << nl;
     }
 
+    Info << "Total omega : " << totalOmega << endl;
+
     Info<< endl;
 
     coeffs_.readIfPresent("useSolarLoad", useSolarLoad_);
 
     if (useSolarLoad_)
     {
+        if (useExternalBeam_)
+        {
+            FatalErrorInFunction
+                << "External beam with fvDOM can not be used "
+                << "with the solar load model"
+                << abort(FatalError);
+        }
         const dictionary& solarDict = this->subDict("solarLoadCoeffs");
         solarLoad_.reset(new solarLoad(solarDict, T_));
 
@@ -296,7 +419,11 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     meshOrientation_
     (
         coeffs_.lookupOrDefault<vector>("meshOrientation", Zero)
-    )
+    ),
+    useExternalBeam_(false),
+    spectralDistribution_(),
+    solarCalculator_(),
+    updateTimeIndex_(0)
 {
     initialise();
 }
@@ -392,7 +519,11 @@ Foam::radiation::fvDOM::fvDOM
     meshOrientation_
     (
         coeffs_.lookupOrDefault<vector>("meshOrientation", Zero)
-    )
+    ),
+    useExternalBeam_(false),
+    spectralDistribution_(),
+    solarCalculator_(),
+    updateTimeIndex_(0)
 {
     initialise();
 }
@@ -435,6 +566,33 @@ void Foam::radiation::fvDOM::calculate()
         solarLoad_->calculate();
     }
 
+    if (useExternalBeam_)
+    {
+        switch (solarCalculator_->sunDirectionModel())
+        {
+            case solarCalculator::mSunDirConstant:
+            {
+                break;
+            }
+            case solarCalculator::mSunDirTracking:
+            {
+                label updateIndex = label
+                (
+                    mesh_.time().value()
+                   /solarCalculator_->sunTrackingUpdateInterval()
+                );
+
+                if (updateIndex > updateTimeIndex_)
+                {
+                    Info << "Updating Sun position..." << endl;
+                    updateTimeIndex_ = updateIndex;
+                    updateRaysDir();
+                }
+                break;
+            }
+        }
+    }
+
     // Set rays convergence false
     List<bool> rayIdConv(nRay_, false);
 
@@ -597,4 +755,10 @@ void Foam::radiation::fvDOM::setRayIdLambdaId
 }
 
 
+const Foam::solarCalculator& Foam::radiation::fvDOM::solarCalc() const
+{
+    return solarCalculator_();
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -32,12 +32,19 @@ Group
 Description
 
     Finite Volume Discrete Ordinates Method. Solves the RTE equation for n
-    directions in a participating media, not including scatter.
+    directions in a participating media, not including scatter and reflective
+    walls.
 
     Available absorption models:
+
     - constantAbsorptionEmission
     - greyMeanAbsoprtionEmission
     - wideBandAbsorptionEmission
+    - multiBandAbsorptionEmission
+
+    This model can handle non-grey participating media using
+    multiBandAbsorptionEmission model. Accordingly the BC for rays should
+    be wideBandDiffussive type
 
 Usage
     \verbatim
@@ -50,6 +57,10 @@ Usage
                                     // iteration
             maxIter     4;          // maximum number of iterations
             meshOrientation    (1 1 1); //Mesh orientation used for 2D and 1D
+
+            useSolarLoad      false;
+            useExternalBeam   true;
+            spectralDistribution (2 1);
         }
 
         solverFreq   1; // Number of flow iterations per radiation iteration
@@ -71,6 +82,16 @@ Usage
       - rays geberated in 3-D using the \c nPhi and \c nTheta entries
       - \c meshOrientation vector is not applicable.
 
+    useSolarLoad calculates the primary and diffusive Sun fluxes on walls in
+    addition to the RTE equations
+
+    useExternalBeam add an external collimated beam to the domain. This option
+    is not available if useSolarLoad is true.
+
+    spectralDistribution is the energy spectral distribution of the collimated
+    external beam.
+
+
 SourceFiles
     fvDOM.C
 
@@ -156,6 +177,18 @@ class fvDOM
         //- Mesh orientation vector
         vector meshOrientation_;
 
+        //- Use external parallel irradiation beam
+        bool useExternalBeam_;
+
+        //- Spectral energy distribution for the external beam
+        scalarList spectralDistribution_;
+
+        //- Solar calculator
+        autoPtr<solarCalculator> solarCalculator_;
+
+        //- Update Sun position index
+        label updateTimeIndex_;
+
 
     // Private Member Functions
 
@@ -213,6 +246,16 @@ public:
                 label& lambdaId
             ) const;
 
+
+            //- Rotate rays according to Sun direction
+            void updateRaysDir();
+
+            //- Rotate rays spheric equator to sunDir
+            void rotateInitialRays(const vector& sunDir);
+
+            //- Align closest ray to sunDir
+            void alignClosestRayToSun(const vector& sunDir);
+
             //- Source term component (for power of T^4)
             virtual tmp<volScalarField> Rp() const;
 
@@ -222,6 +265,9 @@ public:
 
         // Access
 
+            //- Solar calculator
+            const solarCalculator& solarCalc() const;
+
             //- Ray intensity for rayI
             inline const radiativeIntensityRay& IRay(const label rayI) const;
 
@@ -276,6 +322,12 @@ public:
 
             //- Use solar load
             inline bool useSolarLoad() const;
+
+            //- Use external beam
+            inline bool useExternalBeam() const;
+
+            //- Energy spectral distribution for external beam
+            inline const scalarList& spectralDistribution() const;
 };
 
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2008-2011 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -134,4 +134,18 @@ inline bool Foam::radiation::fvDOM::useSolarLoad() const
     return useSolarLoad_;
 }
 
+
+inline bool Foam::radiation::fvDOM::useExternalBeam() const
+{
+    return useExternalBeam_;
+}
+
+
+inline  const Foam::scalarList& Foam::radiation::fvDOM::
+spectralDistribution() const
+{
+    return spectralDistribution_;
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -201,6 +201,12 @@ public:
             //- Return the average vector inside the solid angle
             inline const vector& dAve() const;
 
+             //- Return direction
+            inline vector& d();
+
+            //- Return the average vector inside the solid angle
+            inline vector& dAve();
+
             //- Return the number of bands
             inline scalar nLambda() const;
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H

@@ -82,6 +82,18 @@ inline const Foam::vector& Foam::radiation::radiativeIntensityRay::dAve() const
 }
 
 
+inline Foam::vector& Foam::radiation::radiativeIntensityRay::d()
+{
+    return d_;
+}
+
+
+inline Foam::vector& Foam::radiation::radiativeIntensityRay::dAve()
+{
+    return dAve_;
+}
+
+
 inline Foam::scalar Foam::radiation::radiativeIntensityRay::nLambda() const
 {
     return nLambda_;

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.C

@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiBandZoneAbsorptionEmission.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace radiation
+    {
+        defineTypeNameAndDebug(multiBandZoneAbsorptionEmission, 0);
+
+        addToRunTimeSelectionTable
+        (
+            absorptionEmissionModel,
+            multiBandZoneAbsorptionEmission,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::radiation::multiBandZoneAbsorptionEmission::
+multiBandZoneAbsorptionEmission
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    absorptionEmissionModel(dict, mesh),
+    coeffsDict_(dict.subDict(typeName + "Coeffs")),
+    absCoeffs_(maxBands_),
+    emiCoeffs_(maxBands_),
+    nBands_(0),
+    zoneAbsorptivity_(),
+    zoneEmisivity_(),
+    zoneCells_()
+{
+    coeffsDict_.readEntry("absorptivity", absCoeffs_);
+    coeffsDict_.readEntry("emissivity", emiCoeffs_);
+    nBands_ = absCoeffs_.size();
+
+    const dictionary& zoneDict = coeffsDict_.subDict("zones");
+
+    zoneDict.readEntry("absorptivity", zoneAbsorptivity_);
+    zoneDict.readEntry("emissivity", zoneEmisivity_);
+
+    zoneCells_.setSize(zoneAbsorptivity_.size(), -1);
+
+    label i = 0;
+    forAllConstIters(zoneAbsorptivity_, iter)
+    {
+        label zoneID = mesh.cellZones().findZoneID(iter.key());
+        if (zoneID == -1)
+        {
+            FatalErrorInFunction
+                << "Cannot find cellZone " << iter.key() << endl
+                << "Valid cellZones are " << mesh.cellZones().names()
+                << exit(FatalError);
+        }
+        zoneCells_[i++] = zoneID;
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::radiation::multiBandZoneAbsorptionEmission::
+~multiBandZoneAbsorptionEmission()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::multiBandZoneAbsorptionEmission::aCont
+(
+    const label bandI
+) const
+{
+    tmp<volScalarField> ta
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "a",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("a", dimless/dimLength, absCoeffs_[bandI])
+        )
+    );
+
+    volScalarField& a = ta.ref();
+
+    forAll(zoneCells_, zonei)
+    {
+        const cellZone& cZone = mesh().cellZones()[zoneCells_[zonei]];
+
+        tmp<volScalarField> tzoneAbs(a*0.0);
+        volScalarField& zoneAbs = tzoneAbs.ref();
+
+        const scalarList& abs = zoneAbsorptivity_.find(cZone.name())();
+
+        forAll(cZone, i)
+        {
+            label cellId = cZone[i];
+            zoneAbs[cellId] = abs[bandI] - absCoeffs_[bandI];
+        }
+
+        a += zoneAbs;
+    }
+
+    return ta;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::multiBandZoneAbsorptionEmission::eCont
+(
+    const label bandI
+) const
+{
+    tmp<volScalarField> te
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "e",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("e", dimless/dimLength, emiCoeffs_[bandI])
+        )
+    );
+
+    volScalarField& e = te.ref();
+
+    forAll(zoneCells_, zonei)
+    {
+        const cellZone& cZone = mesh().cellZones()[zoneCells_[zonei]];
+
+        tmp<volScalarField> tzoneEm(e*0.0);
+        volScalarField& zoneEm = tzoneEm.ref();
+
+        const scalarList& emi = zoneEmisivity_.find(cZone.name())();
+
+        forAll(cZone, i)
+        {
+            label cellId = cZone[i];
+            zoneEm[cellId] = emi[bandI] - emiCoeffs_[bandI];
+        }
+        e += zoneEm;
+    }
+
+
+    return te;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::radiation::multiBandZoneAbsorptionEmission::ECont
+(
+    const label bandI
+) const
+{
+    tmp<volScalarField> E
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "E",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
+        )
+    );
+
+    return E;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::radiation::multiBandZoneAbsorptionEmission
+
+Group
+    grpRadiationAbsorptionEmissionSubModels
+
+Description
+
+    multiBandZoneAbsorptionEmission radiation absorption/emission for solids.
+
+
+SourceFiles
+    multiBandZoneAbsorptionEmission.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiBandZoneAbsorptionEmission_H
+#define multiBandZoneAbsorptionEmission_H
+
+#include "absorptionEmissionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class multiBandZoneAbsorptionEmission Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiBandZoneAbsorptionEmission
+:
+    public absorptionEmissionModel
+{
+public:
+
+    // Public data
+
+        //- Maximum number of bands
+        static const int maxBands_ = 5;
+
+
+private:
+
+    // Private data
+
+        //- Absorption model dictionary
+        dictionary coeffsDict_;
+
+        //-  Absorption coefficients outside cellSets
+        scalarList absCoeffs_;
+
+        //-  Emissivity coefficients outside cellSets
+        scalarList emiCoeffs_;
+
+        //- Bands
+        label nBands_;
+
+        //- Cell zones absorptivity
+        HashTable<scalarList> zoneAbsorptivity_;
+
+        //- Cell zones emisivity
+        HashTable<scalarList> zoneEmisivity_;
+
+        //- Cells for each zone
+        labelList zoneCells_;
+
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiBandZoneAbsorptionEmission");
+
+
+    // Constructors
+
+        //- Construct from components
+        multiBandZoneAbsorptionEmission
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~multiBandZoneAbsorptionEmission();
+
+
+    // Member Functions
+
+        // Access
+
+            // Absorption coefficient
+
+                //- Absorption coefficient
+                tmp<volScalarField> aCont(const label bandI) const;
+
+
+            // Emission coefficient
+
+                //- Emission coefficient
+                tmp<volScalarField> eCont(const label bandI) const;
+
+
+            // Emission contribution
+
+                //- Emission contribution
+                tmp<volScalarField> ECont(const label bandI) const;
+
+
+        inline bool isGrey() const
+        {
+            return false;
+        }
+
+        //- Number of bands
+        inline label nBands() const
+        {
+            return nBands_;
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.H

@@ -264,6 +264,12 @@ public:
                 return diffuseSolarRad_;
             }
 
+            //- Return diffuse solar irradiation
+            const scalar& diffuseSolarRad() const
+            {
+                return diffuseSolarRad_;
+            }
+
             //- Return C constant
             scalar C()
             {