View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0001636OpenFOAMBugpublic2015-03-25 14:192015-12-03 15:35
Reportertniemi Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Summary0001636: Questions/issues related to lagrangian radiation models and fvDOM
DescriptionI have several questions/issues related to lagrangian radiation models and fvDOM:

1. Should the radiative heat transfer in ThermoParcel.C be conservative? The reason for asking is that I don't understand the modifications done to the ap and bp terms in calcHeatTransfer. Should the dhsTrans caused by radiation exactly correspond to the radAreaPT4() computations in calc? To me it seems that there might be a problem, as in calc the radiation is computed from the new temperature instead of the old/average temperature used in calcHeatTransfer? I'm using analytic integration.

Anyhow I had some conservation issues, but I was able to eliminate them by assuming a constant parcel temperature for radiation and using T0 temperature in calc when computing radArea terms. I have attached my modifications. (radArea calculation should be modified in every parcel above Thermo) I understand that in case of large temperature changes, this approximation may not be very good, but at least it is conservative.

2. In radiativeIntensityRay.C, correct(), the ECont term is divided with 4 in the IiEq. However, is this division necessary, because if I understand correctly, ECont should already represent 1/4 of the emissions? (in P1.C E_ is multiplied with 4 to get everything)

3. Is it so that the fvDOM radiation model does not work with lagrangian radiation? It seems to only use absorption terms and ignores EDisp. Also Rp and Ru terms do not seem to be correct?

I have attached a modified version, which takes into account the particle emissions in ray equations and also removes parcel absorption/emission from Rp and Ru terms. The modifications seem to work ok in my test cases.
TagsNo tags attached.

Activities

tniemi

2015-03-25 14:19

reporter  

fvDOM.C (14,314 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
     \\/     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 "fvDOM.H"
#include "absorptionEmissionModel.H"
#include "scatterModel.H"
#include "constants.H"
#include "fvm.H"
#include "addToRunTimeSelectionTable.H"

using namespace Foam::constant;
using namespace Foam::constant::mathematical;

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    namespace radiation
    {
        defineTypeNameAndDebug(fvDOM, 0);
        addToRadiationRunTimeSelectionTables(fvDOM);
    }
}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::radiation::fvDOM::initialise()
{
    // 3D
    if (mesh_.nSolutionD() == 3)
    {
        nRay_ = 4*nPhi_*nTheta_;
        IRay_.setSize(nRay_);
        scalar deltaPhi = pi/(2.0*nPhi_);
        scalar deltaTheta = pi/nTheta_;
        label i = 0;
        for (label n = 1; n <= nTheta_; n++)
        {
            for (label m = 1; m <= 4*nPhi_; m++)
            {
                scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
                scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
                IRay_.set
                (
                    i,
                    new radiativeIntensityRay
                    (
                        *this,
                        mesh_,
                        phii,
                        thetai,
                        deltaPhi,
                        deltaTheta,
                        nLambda_,
                        absorptionEmission_,
                        blackBody_,
                        i
                    )
                );
                i++;
            }
        }
    }
    // 2D
    else if (mesh_.nSolutionD() == 2)
    {
        // Currently 2D solution is limited to the x-y plane
        if (mesh_.solutionD()[vector::Z] != -1)
        {
            FatalErrorIn("fvDOM::initialise()")
                << "Currently 2D solution is limited to the x-y plane"
                << exit(FatalError);
        }

        scalar thetai = piByTwo;
        scalar deltaTheta = pi;
        nRay_ = 4*nPhi_;
        IRay_.setSize(nRay_);
        scalar deltaPhi = pi/(2.0*nPhi_);
        label i = 0;
        for (label m = 1; m <= 4*nPhi_; m++)
        {
            scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
            IRay_.set
            (
                i,
                new radiativeIntensityRay
                (
                    *this,
                    mesh_,
                    phii,
                    thetai,
                    deltaPhi,
                    deltaTheta,
                    nLambda_,
                    absorptionEmission_,
                    blackBody_,
                    i
                )
            );
            i++;
        }
    }
    // 1D
    else
    {
        // Currently 1D solution is limited to the x-direction
        if (mesh_.solutionD()[vector::X] != 1)
        {
            FatalErrorIn("fvDOM::initialise()")
                << "Currently 1D solution is limited to the x-direction"
                << exit(FatalError);
        }

        scalar thetai = piByTwo;
        scalar deltaTheta = pi;
        nRay_ = 2;
        IRay_.setSize(nRay_);
        scalar deltaPhi = pi;
        label i = 0;
        for (label m = 1; m <= 2; m++)
        {
            scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
            IRay_.set
            (
                i,
                new radiativeIntensityRay
                (
                    *this,
                    mesh_,
                    phii,
                    thetai,
                    deltaPhi,
                    deltaTheta,
                    nLambda_,
                    absorptionEmission_,
                    blackBody_,
                    i
                )
            );
            i++;
        }
    }


    // Construct absorption field for each wavelength
    forAll(aLambda_, lambdaI)
    {
        aLambda_.set
        (
            lambdaI,
            new volScalarField
            (
                IOobject
                (
                    "aLambda_" + Foam::name(lambdaI) ,
                    mesh_.time().timeName(),
                    mesh_,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                a_
            )
        );
    }

    Info<< "fvDOM : Allocated " << IRay_.size()
        << " rays with average orientation:" << nl;

    if (cacheDiv_)
    {
        Info<< "Caching div fvMatrix..."<< endl;
        for (label lambdaI = 0; lambdaI < nLambda_; lambdaI++)
        {
            fvRayDiv_[lambdaI].setSize(nRay_);

            forAll(IRay_, rayId)
            {
                const surfaceScalarField Ji(IRay_[rayId].dAve() & mesh_.Sf());
                const volScalarField& iRayLambdaI =
                    IRay_[rayId].ILambda(lambdaI);

                fvRayDiv_[lambdaI].set
                (
                    rayId,
                    new fvScalarMatrix
                    (
                        fvm::div(Ji, iRayLambdaI, "div(Ji,Ii_h)")
                    )
                );
            }
        }
    }

    forAll(IRay_, rayId)
    {
        if (omegaMax_ <  IRay_[rayId].omega())
        {
            omegaMax_ = IRay_[rayId].omega();
        }
        Info<< '\t' << IRay_[rayId].I().name() << " : " << "omega : "
            << '\t' << IRay_[rayId].omega() << nl;
    }

    Info<< endl;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
:
    radiationModel(typeName, T),
    G_
    (
        IOobject
        (
            "G",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
    ),
    Qr_
    (
        IOobject
        (
            "Qr",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
    ),
    Qem_
    (
        IOobject
        (
            "Qem",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
    ),
    Qin_
    (
        IOobject
        (
            "Qin",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
    ),
    a_
    (
        IOobject
        (
            "a",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("a", dimless/dimLength, 0.0)
    ),
    nTheta_(readLabel(coeffs_.lookup("nTheta"))),
    nPhi_(readLabel(coeffs_.lookup("nPhi"))),
    nRay_(0),
    nLambda_(absorptionEmission_->nBands()),
    aLambda_(nLambda_),
    blackBody_(nLambda_, T),
    IRay_(0),
    convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
    maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
    fvRayDiv_(nLambda_),
    cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
    omegaMax_(0)
{
    initialise();
}


Foam::radiation::fvDOM::fvDOM
(
    const dictionary& dict,
    const volScalarField& T
)
:
    radiationModel(typeName, dict, T),
    G_
    (
        IOobject
        (
            "G",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
    ),
    Qr_
    (
        IOobject
        (
            "Qr",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
    ),
    Qem_
    (
        IOobject
        (
            "Qem",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
    ),
    Qin_
    (
        IOobject
        (
            "Qin",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
    ),
    a_
    (
        IOobject
        (
            "a",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("a", dimless/dimLength, 0.0)
    ),
    nTheta_(readLabel(coeffs_.lookup("nTheta"))),
    nPhi_(readLabel(coeffs_.lookup("nPhi"))),
    nRay_(0),
    nLambda_(absorptionEmission_->nBands()),
    aLambda_(nLambda_),
    blackBody_(nLambda_, T),
    IRay_(0),
    convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
    maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
    fvRayDiv_(nLambda_),
    cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
    omegaMax_(0)
{
    initialise();
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::radiation::fvDOM::~fvDOM()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

bool Foam::radiation::fvDOM::read()
{
    if (radiationModel::read())
    {
        // Only reading solution parameters - not changing ray geometry
        coeffs_.readIfPresent("convergence", convergence_);
        coeffs_.readIfPresent("maxIter", maxIter_);

        return true;
    }
    else
    {
        return false;
    }
}


void Foam::radiation::fvDOM::calculate()
{
    absorptionEmission_->correct(a_, aLambda_);

    updateBlackBodyEmission();

    // Set rays convergence false
    List<bool> rayIdConv(nRay_, false);

    scalar maxResidual = 0.0;
    label radIter = 0;
    do
    {
        Info<< "Radiation solver iter: " << radIter << endl;

        radIter++;
        maxResidual = 0.0;
        forAll(IRay_, rayI)
        {
            if (!rayIdConv[rayI])
            {
                scalar maxBandResidual = IRay_[rayI].correct();
                maxResidual = max(maxBandResidual, maxResidual);

                if (maxBandResidual < convergence_)
                {
                    rayIdConv[rayI] = true;
                }
            }
        }

    } while (maxResidual > convergence_ && radIter < maxIter_);

    updateG();
}


Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
{
    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "Rp",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            // Only count continuous phase emission
            4.0*absorptionEmission_->aCont()*physicoChemical::sigma 
        )
    );
}


Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::radiation::fvDOM::Ru() const
{

    const DimensionedField<scalar, volMesh>& G =
        G_.dimensionedInternalField();
    const DimensionedField<scalar, volMesh> E =
        absorptionEmission_->ECont()().dimensionedInternalField();
    //const DimensionedField<scalar, volMesh> a =
    //    a_.dimensionedInternalField();
    // Only count continuous phase absorption
    const DimensionedField<scalar, volMesh> a =
    absorptionEmission_->aCont()().dimensionedInternalField();
    return  a*G - E;
}


void Foam::radiation::fvDOM::updateBlackBodyEmission()
{
    for (label j=0; j < nLambda_; j++)
    {
        blackBody_.correct(j, absorptionEmission_->bands(j));
    }
}


void Foam::radiation::fvDOM::updateG()
{
    G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
    Qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
    Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
    Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);

    forAll(IRay_, rayI)
    {
        IRay_[rayI].addIntensity();
        G_ += IRay_[rayI].I()*IRay_[rayI].omega();
        Qr_.boundaryField() += IRay_[rayI].Qr().boundaryField();
        Qem_.boundaryField() += IRay_[rayI].Qem().boundaryField();
        Qin_.boundaryField() += IRay_[rayI].Qin().boundaryField();
    }
}


void Foam::radiation::fvDOM::setRayIdLambdaId
(
    const word& name,
    label& rayId,
    label& lambdaId
) const
{
    // assuming name is in the form: CHARS_rayId_lambdaId
    size_type i1 = name.find_first_of("_");
    size_type i2 = name.find_last_of("_");

    rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
    lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
}


// ************************************************************************* //
fvDOM.C (14,314 bytes)   

tniemi

2015-03-25 14:20

reporter  

radiativeIntensityRay.C (7,717 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
     \\/     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 "radiativeIntensityRay.H"
#include "fvm.H"
#include "fvDOM.H"
#include "constants.H"

using namespace Foam::constant;


const Foam::word
Foam::radiation::radiativeIntensityRay::intensityPrefix("ILambda");


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
(
    const fvDOM& dom,
    const fvMesh& mesh,
    const scalar phi,
    const scalar theta,
    const scalar deltaPhi,
    const scalar deltaTheta,
    const label nLambda,
    const absorptionEmissionModel& absorptionEmission,
    const blackBodyEmission& blackBody,
    const label rayId
)
:
    dom_(dom),
    mesh_(mesh),
    absorptionEmission_(absorptionEmission),
    blackBody_(blackBody),
    I_
    (
        IOobject
        (
            "I" + name(rayId),
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("I", dimMass/pow3(dimTime), 0.0)
    ),
    Qr_
    (
        IOobject
        (
            "Qr" + name(rayId),
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
    ),
    Qin_
    (
        IOobject
        (
            "Qin" + name(rayId),
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
    ),
    Qem_
    (
        IOobject
        (
            "Qem" + name(rayId),
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
    ),
    d_(vector::zero),
    dAve_(vector::zero),
    theta_(theta),
    phi_(phi),
    omega_(0.0),
    nLambda_(nLambda),
    ILambda_(nLambda),
    myRayId_(rayId)
{
    scalar sinTheta = Foam::sin(theta);
    scalar cosTheta = Foam::cos(theta);
    scalar sinPhi = Foam::sin(phi);
    scalar cosPhi = Foam::cos(phi);

    omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
    d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
    dAve_ = vector
    (
        sinPhi
       *Foam::sin(0.5*deltaPhi)
       *(deltaTheta - Foam::cos(2.0*theta)
       *Foam::sin(deltaTheta)),
        cosPhi
       *Foam::sin(0.5*deltaPhi)
       *(deltaTheta - Foam::cos(2.0*theta)
       *Foam::sin(deltaTheta)),
        0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
    );

    autoPtr<volScalarField> IDefaultPtr;

    forAll(ILambda_, lambdaI)
    {
        IOobject IHeader
        (
            intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
            mesh_.time().timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        );

        // check if field exists and can be read
        if (IHeader.headerOk())
        {
            ILambda_.set
            (
                lambdaI,
                new volScalarField(IHeader, mesh_)
            );
        }
        else
        {
            // Demand driven load the IDefault field
            if (!IDefaultPtr.valid())
            {
                IDefaultPtr.reset
                (
                    new volScalarField
                    (
                        IOobject
                        (
                            "IDefault",
                            mesh_.time().timeName(),
                            mesh_,
                            IOobject::MUST_READ,
                            IOobject::NO_WRITE
                        ),
                        mesh_
                    )
                );
            }

            // Reset the MUST_READ flag
            IOobject noReadHeader(IHeader);
            noReadHeader.readOpt() = IOobject::NO_READ;

            ILambda_.set
            (
                lambdaI,
                new volScalarField(noReadHeader, IDefaultPtr())
            );
        }
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::radiation::radiativeIntensityRay::~radiativeIntensityRay()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
{
    // reset boundary heat flux to zero
    Qr_.boundaryField() = 0.0;

    scalar maxResidual = -GREAT;

    forAll(ILambda_, lambdaI)
    {
        const volScalarField& k = dom_.aLambda(lambdaI);

        tmp<fvScalarMatrix> IiEq;

        if (!dom_.cacheDiv())
        {
            const surfaceScalarField Ji(dAve_ & mesh_.Sf());

            IiEq =
            (
                fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
              + fvm::Sp(k*omega_, ILambda_[lambdaI])
            ==
                1.0/constant::mathematical::pi*omega_
              * (
                  // Remove aDisp from k
                    (k-absorptionEmission_.aDisp(lambdaI))*blackBody_.bLambda(lambdaI)
                  + absorptionEmission_.ECont(lambdaI)
                  // Add EDisp term from parcels
                  + absorptionEmission_.EDisp(lambdaI)
                )
            );
        }
        else
        {
            IiEq =
            (
               dom_.fvRayDiv(myRayId_, lambdaI)
             + fvm::Sp(k*omega_, ILambda_[lambdaI])
           ==
               1.0/constant::mathematical::pi*omega_
             * (
                 // Remove aDisp from k
                   (k-absorptionEmission_.aDisp(lambdaI))*blackBody_.bLambda(lambdaI)
                 + absorptionEmission_.ECont(lambdaI)
                 // Add EDisp term from parcels
                 + absorptionEmission_.EDisp(lambdaI)
               )
            );
        }

        IiEq().relax();

        const solverPerformance ILambdaSol = solve
        (
            IiEq(),
            mesh_.solver("Ii")
        );

        const scalar initialRes =
            ILambdaSol.initialResidual()*omega_/dom_.omegaMax();

        maxResidual = max(initialRes, maxResidual);
    }

    return maxResidual;
}


void Foam::radiation::radiativeIntensityRay::addIntensity()
{
    I_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);

    forAll(ILambda_, lambdaI)
    {
        I_ += ILambda_[lambdaI];
    }
}


// ************************************************************************* //
radiativeIntensityRay.C (7,717 bytes)   

tniemi

2015-03-25 14:20

reporter  

ThermoParcel.C (10,001 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     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 "ThermoParcel.H"
#include "physicoChemicalConstants.H"

using namespace Foam::constant;

// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //

template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::setCellValues
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    ParcelType::setCellValues(td, dt, cellI);

    tetIndices tetIs = this->currentTetIndices();

    Cpc_ = td.CpInterp().interpolate(this->position(), tetIs);

    Tc_ = td.TInterp().interpolate(this->position(), tetIs);

    if (Tc_ < td.cloud().constProps().TMin())
    {
        if (debug)
        {
            WarningIn
            (
                "void Foam::ThermoParcel<ParcelType>::setCellValues"
                "("
                    "TrackData&, "
                    "const scalar, "
                    "const label"
                ")"
            )   << "Limiting observed temperature in cell " << cellI << " to "
                << td.cloud().constProps().TMin() <<  nl << endl;
        }

        Tc_ = td.cloud().constProps().TMin();
    }
}


template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    this->Uc_ += td.cloud().UTrans()[cellI]/this->massCell(cellI);

    const scalar CpMean = td.CpInterp().psi()[cellI];
    Tc_ += td.cloud().hsTrans()[cellI]/(CpMean*this->massCell(cellI));

    if (Tc_ < td.cloud().constProps().TMin())
    {
        if (debug)
        {
            WarningIn
            (
                "void Foam::ThermoParcel<ParcelType>::cellValueSourceCorrection"
                "("
                    "TrackData&, "
                    "const scalar, "
                    "const label"
                ")"
            )   << "Limiting observed temperature in cell " << cellI << " to "
                << td.cloud().constProps().TMin() <<  nl << endl;
        }

        Tc_ = td.cloud().constProps().TMin();
    }
}


template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
(
    TrackData& td,
    const label cellI,
    const scalar T,
    scalar& Ts,
    scalar& rhos,
    scalar& mus,
    scalar& Pr,
    scalar& kappas
) const
{
    // Surface temperature using two thirds rule
    Ts = (2.0*T + Tc_)/3.0;

    if (Ts < td.cloud().constProps().TMin())
    {
        if (debug)
        {
            WarningIn
            (
                "void Foam::ThermoParcel<ParcelType>::calcSurfaceValues"
                "("
                    "TrackData&, "
                    "const label, "
                    "const scalar, "
                    "scalar&, "
                    "scalar&, "
                    "scalar&, "
                    "scalar&, "
                    "scalar&"
                ") const"
            )   << "Limiting parcel surface temperature to "
                << td.cloud().constProps().TMin() <<  nl << endl;
        }

        Ts = td.cloud().constProps().TMin();
    }

    // Assuming thermo props vary linearly with T for small d(T)
    const scalar TRatio = Tc_/Ts;

    rhos = this->rhoc_*TRatio;

    tetIndices tetIs = this->currentTetIndices();
    mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
    kappas = td.kappaInterp().interpolate(this->position(), tetIs)/TRatio;

    Pr = Cpc_*mus/kappas;
    Pr = max(ROOTVSMALL, Pr);
}


template<class ParcelType>
template<class TrackData>
void Foam::ThermoParcel<ParcelType>::calc
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    // Define local properties at beginning of time step
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    const scalar np0 = this->nParticle_;
    const scalar mass0 = this->mass();


    // Calc surface values
    // ~~~~~~~~~~~~~~~~~~~
    scalar Ts, rhos, mus, Pr, kappas;
    calcSurfaceValues(td, cellI, this->T_, Ts, rhos, mus, Pr, kappas);

    // Reynolds number
    scalar Re = this->Re(this->U_, this->d_, rhos, mus);


    // Sources
    // ~~~~~~~

    // Explicit momentum source for particle
    vector Su = vector::zero;

    // Linearised momentum source coefficient
    scalar Spu = 0.0;

    // Momentum transfer from the particle to the carrier phase
    vector dUTrans = vector::zero;

    // Explicit enthalpy source for particle
    scalar Sh = 0.0;

    // Linearised enthalpy source coefficient
    scalar Sph = 0.0;

    // Sensible enthalpy transfer from the particle to the carrier phase
    scalar dhsTrans = 0.0;


    // Heat transfer
    // ~~~~~~~~~~~~~

    // Sum Ni*Cpi*Wi of emission species
    scalar NCpW = 0.0;

    const scalar T0 = this->T_;
    // Calculate new particle temperature
    this->T_ =
        this->calcHeatTransfer
        (
            td,
            dt,
            cellI,
            Re,
            Pr,
            kappas,
            NCpW,
            Sh,
            dhsTrans,
            Sph
        );


    // Motion
    // ~~~~~~

    // Calculate new particle velocity
    this->U_ =
        this->calcVelocity(td, dt, cellI, Re, mus, mass0, Su, dUTrans, Spu);


    //  Accumulate carrier phase source terms
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if (td.cloud().solution().coupled())
    {
        // Update momentum transfer
        td.cloud().UTrans()[cellI] += np0*dUTrans;

        // Update momentum transfer coefficient
        td.cloud().UCoeff()[cellI] += np0*Spu;

        // Update sensible enthalpy transfer
        td.cloud().hsTrans()[cellI] += np0*dhsTrans;

        // Update sensible enthalpy coefficient
        td.cloud().hsCoeff()[cellI] += np0*Sph;

        // Update radiation fields
        if (td.cloud().radiation())
        {
            const scalar ap = this->areaP();
            //Use T0
            const scalar T4 = pow4(T0);
            td.cloud().radAreaP()[cellI] += dt*np0*ap;
            td.cloud().radT4()[cellI] += dt*np0*T4;
            td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
        }
    }
}


template<class ParcelType>
template<class TrackData>
Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
(
    TrackData& td,
    const scalar dt,
    const label cellI,
    const scalar Re,
    const scalar Pr,
    const scalar kappa,
    const scalar NCpW,
    const scalar Sh,
    scalar& dhsTrans,
    scalar& Sph
)
{
    if (!td.cloud().heatTransfer().active())
    {
        return T_;
    }

    const scalar d = this->d();
    const scalar rho = this->rho();

    // Calc heat transfer coefficient
    scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);

    if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
    {
        return
            max
            (
                T_ + dt*Sh/(this->volume(d)*rho*Cp_),
                td.cloud().constProps().TMin()
            );
    }

    htc = max(htc, ROOTVSMALL);
    const scalar As = this->areaS(d);

    scalar ap = Tc_ + Sh/As/htc;
    scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
    if (td.cloud().radiation())
    {
        tetIndices tetIs = this->currentTetIndices();
        const scalar Gc = td.GInterp().interpolate(this->position(), tetIs);
        const scalar sigma = physicoChemical::sigma.value();
        const scalar epsilon = td.cloud().constProps().epsilon0();

        //ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T_)/htc);
        //Assume constant source epsilon*sigma*(Gc/(4.0)-pow4(T_))
        ap = (ap + epsilon/htc*(Gc/(4.0)-sigma*pow4(T_)));
        bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T_)));
    }
    bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;

    // Integrate to find the new parcel temperature
    IntegrationScheme<scalar>::integrationResult Tres =
        td.cloud().TIntegrator().integrate(T_, dt, ap*bp, bp);

    scalar Tnew =
        min
        (
            max
            (
                Tres.value(),
                td.cloud().constProps().TMin()
            ),
            td.cloud().constProps().TMax()
        );

    Sph = dt*htc*As;

    dhsTrans += Sph*(Tres.average() - Tc_);

    return Tnew;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class ParcelType>
Foam::ThermoParcel<ParcelType>::ThermoParcel
(
    const ThermoParcel<ParcelType>& p
)
:
    ParcelType(p),
    T_(p.T_),
    Cp_(p.Cp_),
    Tc_(p.Tc_),
    Cpc_(p.Cpc_)
{}


template<class ParcelType>
Foam::ThermoParcel<ParcelType>::ThermoParcel
(
    const ThermoParcel<ParcelType>& p,
    const polyMesh& mesh
)
:
    ParcelType(p, mesh),
    T_(p.T_),
    Cp_(p.Cp_),
    Tc_(p.Tc_),
    Cpc_(p.Cpc_)
{}


// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //

#include "ThermoParcelIO.C"

// ************************************************************************* //
ThermoParcel.C (10,001 bytes)   

henry

2015-03-27 15:25

manager   ~0004498

Lagrangian has not been tested with fvDOM and there are no tutorial examples.

Thanks for the patches, I will study them and get back to you. It would be very useful to have a simple tutorial case to test radiative heat transfer with Lagrangian tracking; could you propose or provide something suitable?

henry

2015-05-26 11:48

manager   ~0004812

I am looking into including these patches but ideally need a test-case, do you have something suitable?

henry

2015-06-12 15:43

manager   ~0004917

> 1. Should the radiative heat transfer in ThermoParcel.C be conservative?

Having conservation at every iteration would be good unless it compromises convergence. Often a non-conservative semi-implicit approach is applied for stability and conservation is achieved at convergence which may require several iterations. If there is no stability issue then a conservative explicit approach is fine.

> 2. In radiativeIntensityRay.C

I think your analysis is correct.

> 3. Is it so that the fvDOM radiation model does not work with lagrangian
> radiation?

Correct, I have applied your changes to OpenFOAM-dev:
    commit d8f45cf8c1232f9be0932501148c1416e77011ef

I do not have a test-case to check the merge of these patches, could you test and report back?

tniemi

2015-06-15 14:11

reporter   ~0004925

Hi

Sorry for my inactivity, I have been very busy doing other stuff.

> I do not have a test-case to check the merge of these patches, could you test and report back?

I have checked the patches and they seem to be ok. Unfortunately, I don't have a good stand-alone test case either, but we have used this formulation for a while without any issues. So in my opinion the commit is ok.

However, the modification "this->T_ -> T0" in the "Update radiation fields" step in ThermoParcel.C calc should also be applied to ReactingParcel and ReactingMultiphaseParcel, since they also do the same calculation in their calc funtions.

henry

2015-06-15 17:43

manager   ~0004934

Resolved by commit 6f5936b373a5e5456eae2f5ad79d4472e95eee2b

karlvirgil

2015-11-30 15:03

reporter   ~0005691

In radiativeIntensityRay.C, if the absorptionEmission_.ECont is not divided by 4, then the incident radiation flux on a surface, Qin, is 4 times larger than previously. This goes against all of the validation that we have done over the last 8 years for flame radiation on surfaces.

+ absorptionEmission_.ECont(lambdaI)/4

tniemi

2015-12-01 09:07

reporter   ~0005698

Last edited: 2015-12-01 09:11

 Hi karlvirgil,

In order to get the same results as before, you would have to redefine your ECont to be 1/4 th of the value it was before the change, as ECont is now defined to be 1/4 th of the total emission.

The reason why I suggested this change was that it would make the P1 and fvDOM models consistent. In P1 model the E is multiplied by four and it is therefore assumed that E is defined as 1/4 th of the total emissions. However, in the original fvDOM implementation E was assumed to represent the total emissions and for that reason it was divided by four (the division is required, because the sum of omegas over all rays is 4*pi). This is ok, but the different definition of E meant that it would not be possible to use the same ECont model for both radiation models.

One could argue, that the it would be more logical to define ECont without the 1/4 th division. However, then the P1 model would have to be changed. Also the lagrangian cloudAbsorptionEmission currently assumes that EDisp represents 1/4 th of the total emissions.

However, as I looked through the code I noticed a bug that the change has caused, which I didn't notice before. In fvDOM.C in Ru(), the expression should now be a*G - 4.0*E; (the same as in P1), because E is only 1/4 th of the total emissions now.

tniemi

2015-12-01 14:46

reporter   ~0005705

Hi

karlvirgil, in your validation cases, have you used either greyMeanSolidAbsorptionEmission or wideBandAbsorptionEmission models or have you used your own custom absorbtionEmission model? I realized that the greyMeanSolidAbsorptionEmission and wideBandAbsorptionEmission models can use ECont and those models also assume that ECont represents the total emission per surface area. This means, that they assume that ECont should be divided by 4.

To summarize I think there are now two problems: the greyMeanSolidAbsorptionEmission and wideBandAbsorptionEmission use a different definition of ECont than the radiation models and the ECont term in fvDOM Ru function is currently missing the multiplier.

In order to fix these, I can come up with 3 alternatives:

1. Keep ECont defined as it is now, meaning it is emission per projected area/cross sectional area. Fix the bug in fvDOM Ru (multiply ECont by 4) and modify the greyMeanSolidAbsorptionEmission and wideBandAbsorptionEmission so that ECont is divided by 4 there. No modification necessary for P1 or lagrangian models.

2. Define ECont to be the total emission per surface area as it was before in fvDOM. Divide it by 4 in fvDOM, split the E term in P1 and move out the ECont part so that it is not multiplied by 4 (EDisp has to be multiplied). Remove also the multiplication from the Ru function in P1.

3. Define both ECont and EDisp to be the total emission per surface area, apply multiplication by 4 in cloudAbsorptionEmission model (the only place that uses EDisp?). Do not multiply E in P1 at all, divide both ECont and EDisp in fvDOM.

Personally for me it doesn't matter how the ECont is defined, but in my opinion it would be good that ECont has the same definition in both radiation models. Also in options 1 and 3 there is the additional benefit that ECont and EDisp have the same units and can be divided/multiplied in the same way.

henry

2015-12-01 16:24

manager   ~0005708

I have applied option 3 provided by Timo to OpenFOAM-dev:

commit a3301464bd6b110b8e32fa8e123df4c55548d420

Let me know if this is satisfartory and if so I will also apply the correction to OpenFOAM-3.0.x.

karlvirgil

2015-12-01 18:24

reporter   ~0005710

Hi tniemi and henry.

Ankur, my colleague who is well versed in the radiation portions of OpenFOAM, will be responding to this thread shortly with what we perceive as the best solution.

user1287

2015-12-01 19:53

  ~0005711

Hello,

It would be logical to define Econt and EDisp as the total emission per unit volume, as suggested in the option 3 in tniemi's post. The definition of ECont and EDisp as total emission divided by 4 is rather confusing, and would make the code hard to follow.

-Ankur

henry

2015-12-01 22:18

manager   ~0005712

I have already applied option 3 provided by Timo to OpenFOAM-dev, could you please test it and let us know if it is satisfartory?

karlvirgil

2015-12-02 13:50

reporter   ~0005717

I'll test it today and let you know.

karlvirgil

2015-12-03 14:40

reporter   ~0005719

After running a few tests, these changes appear to be satisfactory.

henry

2015-12-03 15:35

manager   ~0005720

Resolved in OpenFOAM-3.0.x by commit df1ea7d5e5035a2628e20ef4b08527f9d4f9a8f7
Resolved in OpenFOAM-dev by commit a3301464bd6b110b8e32fa8e123df4c55548d420

Issue History

Date Modified Username Field Change
2015-03-25 14:19 tniemi New Issue
2015-03-25 14:19 tniemi File Added: fvDOM.C
2015-03-25 14:20 tniemi File Added: radiativeIntensityRay.C
2015-03-25 14:20 tniemi File Added: ThermoParcel.C
2015-03-27 15:25 henry Note Added: 0004498
2015-05-26 11:48 henry Note Added: 0004812
2015-06-12 15:43 henry Note Added: 0004917
2015-06-15 14:11 tniemi Note Added: 0004925
2015-06-15 17:43 henry Note Added: 0004934
2015-06-15 17:43 henry Status new => resolved
2015-06-15 17:43 henry Resolution open => fixed
2015-06-15 17:43 henry Assigned To => henry
2015-11-19 21:57 henry Status resolved => feedback
2015-11-19 21:57 henry Resolution fixed => reopened
2015-11-30 15:03 karlvirgil Note Added: 0005691
2015-12-01 09:07 tniemi Note Added: 0005698
2015-12-01 09:07 tniemi Status feedback => assigned
2015-12-01 09:11 tniemi Note Edited: 0005698
2015-12-01 14:46 tniemi Note Added: 0005705
2015-12-01 16:24 henry Note Added: 0005708
2015-12-01 18:24 karlvirgil Note Added: 0005710
2015-12-01 19:53 user1287 Note Added: 0005711
2015-12-01 22:18 henry Note Added: 0005712
2015-12-02 13:50 karlvirgil Note Added: 0005717
2015-12-03 14:40 karlvirgil Note Added: 0005719
2015-12-03 15:35 henry Note Added: 0005720
2015-12-03 15:35 henry Status assigned => resolved
2015-12-03 15:35 henry Resolution reopened => fixed