diff --git a/src/thermophysicalModels/radiation/Make/files b/src/thermophysicalModels/radiation/Make/files index 6c5523bbe..f76487c47 100644 --- a/src/thermophysicalModels/radiation/Make/files +++ b/src/thermophysicalModels/radiation/Make/files @@ -33,6 +33,7 @@ derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C +derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.C derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C LIB = $(FOAM_LIBBIN)/libradiation diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C index b7b13b408..f458c5190 100644 --- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C @@ -29,13 +29,64 @@ License #include "volFields.H" #include "fvDOM.H" -#include "radiationConstants.H" -#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace radiation +{ +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void +greyDiffusiveRadiationMixedFvPatchScalarField::calcSumOutgoingAngles() const +{ + if (sumOutgoingAnglesPtr_) + { + FatalErrorIn + ( + "wideBandDiffusiveRadiationMixedFvPatchScalarField" + "::calcSumOutgoingAngles()" + ) << "sumOutgoingAnglesPtr_ already calculated" + << abort(FatalError); + } + + sumOutgoingAnglesPtr_ = new scalarField(this->size(), 0.); + scalarField& sumOutgoingAngles = *sumOutgoingAnglesPtr_; + + + // Get access to radiation model, and recast as fvDOM + const radiationModel& radiation = + db().lookupObject("radiationProperties"); + + const fvDOM& dom = dynamic_cast(radiation); + + + for(label rayI = 0; rayI < dom.nRay(); rayI++) + { + // Calculate cosine of angle between face and ray + scalarField outgoingAngles = this->patch().nf() & dom.IRay(rayI).dAve(); + + // For outgoing rays, outgoingAngles will be negative + sumOutgoingAngles += neg(outgoingAngles)*(-outgoingAngles); + } +} + +const Foam::scalarField& +greyDiffusiveRadiationMixedFvPatchScalarField::sumOutgoingAngles() const +{ + if (!sumOutgoingAnglesPtr_) + { + calcSumOutgoingAngles(); + } + + return *sumOutgoingAnglesPtr_; +} // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField:: +greyDiffusiveRadiationMixedFvPatchScalarField:: greyDiffusiveRadiationMixedFvPatchScalarField ( const fvPatch& p, @@ -44,7 +95,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(p, iF), TName_("undefinedT"), - emissivity_(0.0) + emissivity_(0.0), + sumOutgoingAnglesPtr_(NULL) { refValue() = 0.0; refGrad() = 0.0; @@ -52,7 +104,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField } -Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField:: +greyDiffusiveRadiationMixedFvPatchScalarField:: greyDiffusiveRadiationMixedFvPatchScalarField ( const greyDiffusiveRadiationMixedFvPatchScalarField& ptf, @@ -63,11 +115,12 @@ greyDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(ptf, p, iF, mapper), TName_(ptf.TName_), - emissivity_(ptf.emissivity_) + emissivity_(ptf.emissivity_), + sumOutgoingAnglesPtr_(NULL) {} -Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField:: +greyDiffusiveRadiationMixedFvPatchScalarField:: greyDiffusiveRadiationMixedFvPatchScalarField ( const fvPatch& p, @@ -77,7 +130,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(p, iF), TName_(dict.lookup("T")), - emissivity_(readScalar(dict.lookup("emissivity"))) + emissivity_(readScalar(dict.lookup("emissivity"))), + sumOutgoingAnglesPtr_(NULL) { if (dict.found("refValue")) { @@ -110,7 +164,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField } -Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField:: +greyDiffusiveRadiationMixedFvPatchScalarField:: greyDiffusiveRadiationMixedFvPatchScalarField ( const greyDiffusiveRadiationMixedFvPatchScalarField& ptf @@ -118,11 +172,12 @@ greyDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(ptf), TName_(ptf.TName_), - emissivity_(ptf.emissivity_) + emissivity_(ptf.emissivity_), + sumOutgoingAnglesPtr_(NULL) {} -Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField:: +greyDiffusiveRadiationMixedFvPatchScalarField:: greyDiffusiveRadiationMixedFvPatchScalarField ( const greyDiffusiveRadiationMixedFvPatchScalarField& ptf, @@ -131,100 +186,88 @@ greyDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(ptf, iF), TName_(ptf.TName_), - emissivity_(ptf.emissivity_) + emissivity_(ptf.emissivity_), + sumOutgoingAnglesPtr_(NULL) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField:: -updateCoeffs() +void greyDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs() { if (this->updated()) { return; } - const scalarField& Tp = - lookupPatchField(TName_); + + + const label patchI = this->patch().index(); + + // Access radiation model const radiationModel& radiation = db().lookupObject("radiationProperties"); - const fvDOM& dom(refCast(radiation)); + const fvDOM& dom = dynamic_cast(radiation); + if (dom.nLambda() == 0) + { + FatalErrorIn + ( + "" + "wideBandDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs" + ) << " a non-grey boundary condition is used with a grey " + << "absorption model" << nl << exit(FatalError); + } + + // Get rayId and lambda Id for this ray label rayId = -1; label lambdaId = -1; dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId); - const label patchI = patch().index(); + // Make shortcut to ray belonging to this field + const radiativeIntensityRay& ray = dom.IRay(rayId); - if (dom.nLambda() != 1) - { - FatalErrorIn + // Access incoming radiation for this patch for this band + const scalarField& Qin = dom.Qin(lambdaId)[patchI]; + + // Access black body radiation for this patch for this band + const scalarField& Eb = + dom.blackBody().bLambda(lambdaId).boundaryField()[patchI]; + + // Get face normals + vectorField nHat = patch().nf(); + + // Calculate cos of incoming angle of current ray with every face + scalarField incomingAngle = this->patch().nf() & ray.dAve(); + + // Set to zeroGradient (=0; incomingAngle > 0) for faces with incoming rays + // and to fixedValue (=1; incomingAngle < 0) for outgoing rays + this->valueFraction() = neg(incomingAngle); + + // Set intensity value for fixedValue part (see reference in header file) + this->refValue() = + emissivity_*Eb + ((1 - emissivity_)*Qin)/sumOutgoingAngles(); + + // Update boundary field now, so values for incoming and outgoing rays + // are in balance + scalarField::operator= + ( + this->valueFraction()*this->refValue() + + + (1.0 - this->valueFraction())* ( - "Foam::radiation::" - "greyDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs" - ) << " a grey boundary condition is used with a non-grey " - << "absorption model" << nl << exit(FatalError); - } - - scalarField& Iw = *this; - vectorField n = patch().Sf()/patch().magSf(); - - radiativeIntensityRay& ray = - const_cast(dom.IRay(rayId)); - - ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve()); - - forAll(Iw, faceI) - { - scalar Ir = 0.0; - - for (label rayI=0; rayI < dom.nRay(); rayI++) - { - const vector& d = dom.IRay(rayI).d(); - - const scalarField& IFace = - dom.IRay(rayI).ILambda(lambdaId).boundaryField()[patchI]; - - if ((-n[faceI] & d) < 0.0) - { - // q into the wall - const vector& dAve = dom.IRay(rayI).dAve(); - Ir += IFace[faceI]*mag(n[faceI] & dAve); - } - } - - const vector& d = dom.IRay(rayId).d(); - - if ((-n[faceI] & d) > 0.0) - { - // direction out of the wall - refGrad()[faceI] = 0.0; - valueFraction()[faceI] = 1.0; - refValue()[faceI] = - ( - Ir*(1.0 - emissivity_) - + emissivity_*radiation::sigmaSB.value()*pow4(Tp[faceI]) - ) - /mathematicalConstant::pi; - - } - else - { - // direction into the wall - valueFraction()[faceI] = 0.0; - refGrad()[faceI] = 0.0; - refValue()[faceI] = 0.0; //not used - } - } + this->patchInternalField() + + this->refGrad()/this->patch().deltaCoeffs() + ) + ); mixedFvPatchScalarField::updateCoeffs(); } -void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::write +void greyDiffusiveRadiationMixedFvPatchScalarField::write ( Ostream& os ) const @@ -237,17 +280,13 @@ void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::write // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -namespace Foam -{ -namespace radiation -{ - makePatchTypeField - ( - fvPatchScalarField, - greyDiffusiveRadiationMixedFvPatchScalarField - ); -} -} +makePatchTypeField +( + fvPatchScalarField, + greyDiffusiveRadiationMixedFvPatchScalarField +); +} // End namespace radiation +} // End namespace Foam // ************************************************************************* // diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H index 515514117..74dd670ab 100644 --- a/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H @@ -25,7 +25,9 @@ Class Foam::greyDiffusiveRadiationMixedFvPatchScalarField Description - Radiation temperature specified + Radiation boundary will diffusive reflection and black body emission. + See Modest "Radiative heat transfer", Capter 16.6 "The finite volume method" + Eq.16.63 SourceFiles greyDiffusiveRadiationMixedFvPatchScalarField.C @@ -59,6 +61,15 @@ class greyDiffusiveRadiationMixedFvPatchScalarField //- Emissivity scalar emissivity_; + //- Sum of cosine factors of all outgoing rays (lazy evaluated) + mutable scalarField *sumOutgoingAnglesPtr_; + + + // Private member functions + + //- Calculate sum of cosine factors of all outgoing rays + void calcSumOutgoingAngles() const; + public: @@ -157,6 +168,9 @@ public: return emissivity_; } + //- Return sum of cosine factors of all outgoing rays + const scalarField& sumOutgoingAngles() const; + // Evaluation functions diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C index 53b749cb1..2716a82be 100644 --- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C @@ -29,13 +29,64 @@ License #include "volFields.H" #include "fvDOM.H" -#include "wideBandAbsorptionEmission.H" -#include "radiationConstants.H" -#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace radiation +{ +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void +wideBandDiffusiveRadiationMixedFvPatchScalarField::calcSumOutgoingAngles() const +{ + if (sumOutgoingAnglesPtr_) + { + FatalErrorIn + ( + "wideBandDiffusiveRadiationMixedFvPatchScalarField" + "::calcSumOutgoingAngles()" + ) << "sumOutgoingAnglesPtr_ already calculated" + << abort(FatalError); + } + + sumOutgoingAnglesPtr_ = new scalarField(this->size(), 0.); + scalarField& sumOutgoingAngles = *sumOutgoingAnglesPtr_; + + + // Get access to radiation model, and recast as fvDOM + const radiationModel& radiation = + db().lookupObject("radiationProperties"); + + const fvDOM& dom = dynamic_cast(radiation); + + + for(label rayI = 0; rayI < dom.nRay(); rayI++) + { + // Calculate cosine of angle between face and ray + scalarField outgoingAngles = this->patch().nf() & dom.IRay(rayI).dAve(); + + // For outgoing rays, outgoingAngles will be negative + sumOutgoingAngles += neg(outgoingAngles)*(-outgoingAngles); + } +} + +const Foam::scalarField& +wideBandDiffusiveRadiationMixedFvPatchScalarField::sumOutgoingAngles() const +{ + if (!sumOutgoingAnglesPtr_) + { + calcSumOutgoingAngles(); + } + + return *sumOutgoingAnglesPtr_; +} + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField:: +wideBandDiffusiveRadiationMixedFvPatchScalarField:: wideBandDiffusiveRadiationMixedFvPatchScalarField ( const fvPatch& p, @@ -44,7 +95,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(p, iF), TName_("undefinedT"), - emissivity_(0.0) + emissivity_(0.0), + sumOutgoingAnglesPtr_(NULL) { refValue() = 0.0; refGrad() = 0.0; @@ -52,7 +104,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField } -Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField:: +wideBandDiffusiveRadiationMixedFvPatchScalarField:: wideBandDiffusiveRadiationMixedFvPatchScalarField ( const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf, @@ -63,11 +115,12 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(ptf, p, iF, mapper), TName_(ptf.TName_), - emissivity_(ptf.emissivity_) + emissivity_(ptf.emissivity_), + sumOutgoingAnglesPtr_(NULL) {} -Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField:: +wideBandDiffusiveRadiationMixedFvPatchScalarField:: wideBandDiffusiveRadiationMixedFvPatchScalarField ( const fvPatch& p, @@ -77,7 +130,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(p, iF), TName_(dict.lookup("T")), - emissivity_(readScalar(dict.lookup("emissivity"))) + emissivity_(readScalar(dict.lookup("emissivity"))), + sumOutgoingAnglesPtr_(NULL) { if (dict.found("value")) { @@ -95,7 +149,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField { // No value given. Restart as fixedValue b.c. - // Bugfix: Do not initialize from temperautre because it is unavailable + // Bugfix: Do not initialize from temperature because it is unavailable // when running, e.g. decomposePar and loading radiation as // shared library. Initialize to zero instead. // 26 Mar 2014 - DC @@ -109,7 +163,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField } -Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField:: +wideBandDiffusiveRadiationMixedFvPatchScalarField:: wideBandDiffusiveRadiationMixedFvPatchScalarField ( const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf @@ -117,11 +171,12 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(ptf), TName_(ptf.TName_), - emissivity_(ptf.emissivity_) + emissivity_(ptf.emissivity_), + sumOutgoingAnglesPtr_(NULL) {} -Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField:: +wideBandDiffusiveRadiationMixedFvPatchScalarField:: wideBandDiffusiveRadiationMixedFvPatchScalarField ( const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf, @@ -130,97 +185,85 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField : mixedFvPatchScalarField(ptf, iF), TName_(ptf.TName_), - emissivity_(ptf.emissivity_) + emissivity_(ptf.emissivity_), + sumOutgoingAnglesPtr_(NULL) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField:: -updateCoeffs() +void wideBandDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs() { if (this->updated()) { return; } + const label patchI = this->patch().index(); + + // Access radiation model const radiationModel& radiation = db().lookupObject("radiationProperties"); - const fvDOM& dom(refCast(radiation)); - - label rayId = -1; - label lambdaId = -1; - dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId); - - const label patchI = patch().index(); + const fvDOM& dom = dynamic_cast(radiation); if (dom.nLambda() == 0) { FatalErrorIn ( - "Foam::radiation::" + "" "wideBandDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs" ) << " a non-grey boundary condition is used with a grey " << "absorption model" << nl << exit(FatalError); } - scalarField& Iw = *this; - vectorField n = patch().Sf()/patch().magSf(); + // Get rayId and lambda Id for this ray + label rayId = -1; + label lambdaId = -1; + dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId); - radiativeIntensityRay& ray = - const_cast(dom.IRay(rayId)); + // Make shortcut to ray belonging to this field + const radiativeIntensityRay& ray = dom.IRay(rayId); - ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve()); + // Access incoming radiation for this patch for this band + const scalarField& Qin = dom.Qin(lambdaId)[patchI]; - const scalarField Eb = + // Access black body radiation for this patch for this band + const scalarField& Eb = dom.blackBody().bLambda(lambdaId).boundaryField()[patchI]; - forAll(Iw, faceI) - { - scalar Ir = 0.0; - for (label rayI=0; rayI < dom.nRay(); rayI++) - { - const vector& d = dom.IRay(rayI).d(); + // Get face normals + vectorField nHat = patch().nf(); - const scalarField& IFace = - dom.IRay(rayI).ILambda(lambdaId).boundaryField()[patchI]; + // Calculate cos of incoming angle of current ray with every face + scalarField incomingAngle = this->patch().nf() & ray.dAve(); - if ((-n[faceI] & d) < 0.0) // qin into the wall - { - const vector& dAve = dom.IRay(rayI).dAve(); - Ir = Ir + IFace[faceI]*mag(n[faceI] & dAve); - } - } + // Set to zeroGradient (=0; incomingAngle > 0) for faces with incoming rays + // and to fixedValue (=1; incomingAngle < 0) for outgoing rays + this->valueFraction() = neg(incomingAngle); - const vector& d = dom.IRay(rayId).d(); + // Set intensity value for fixedValue part (see reference in header file) + this->refValue() = + emissivity_*Eb + ((1 - emissivity_)*Qin)/sumOutgoingAngles(); - if ((-n[faceI] & d) > 0.0) - { - // direction out of the wall - refGrad()[faceI] = 0.0; - valueFraction()[faceI] = 1.0; - refValue()[faceI] = - ( - Ir*(1.0 - emissivity_) - + emissivity_*Eb[faceI] - ) - /mathematicalConstant::pi; - } - else - { - // direction into the wall - valueFraction()[faceI] = 0.0; - refGrad()[faceI] = 0.0; - refValue()[faceI] = 0.0; //not used - } - } + // Update boundary field now, so values for incoming and outgoing rays + // are in balance + scalarField::operator= + ( + this->valueFraction()*this->refValue() + + + (1.0 - this->valueFraction())* + ( + this->patchInternalField() + + this->refGrad()/this->patch().deltaCoeffs() + ) + ); mixedFvPatchScalarField::updateCoeffs(); } -void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::write +void wideBandDiffusiveRadiationMixedFvPatchScalarField::write ( Ostream& os ) const @@ -234,17 +277,14 @@ void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::write // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -namespace Foam -{ -namespace radiation -{ - makePatchTypeField - ( - fvPatchScalarField, - wideBandDiffusiveRadiationMixedFvPatchScalarField - ); -} -} +makePatchTypeField +( + fvPatchScalarField, + wideBandDiffusiveRadiationMixedFvPatchScalarField +); +} // End namespace radiation +} // End namespace Foam + // ************************************************************************* // diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H index 70ced563f..22fd92c4b 100644 --- a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H @@ -25,7 +25,9 @@ Class Foam::wideBandDiffusiveRadiationMixedFvPatchScalarField Description - Radiation temperature specified + Radiation boundary will diffusive reflection and black body emission. + See Modest "Radiative heat transfer", Capter 16.6 "The finite volume method" + Eq.16.63 SourceFiles wideBandDiffusiveRadiationMixedFvPatchScalarField.C @@ -59,6 +61,15 @@ class wideBandDiffusiveRadiationMixedFvPatchScalarField //- Emissivity scalar emissivity_; + //- Sum of cosine factors of all outgoing rays (lazy evaluated) + mutable scalarField *sumOutgoingAnglesPtr_; + + + // Private member functions + + //- Calculate sum of cosine factors of all outgoing rays + void calcSumOutgoingAngles() const; + public: @@ -157,6 +168,9 @@ public: return emissivity_; } + //- Return sum of cosine factors of all outgoing rays + const scalarField& sumOutgoingAngles() const; + // Evaluation functions diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.C new file mode 100644 index 000000000..ef8f858ba --- /dev/null +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.C @@ -0,0 +1,416 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "radiationModel.H" +#include "fvDOM.H" +#include "wideBandSpecularRadiationMixedFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace radiation +{ +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void wideBandSpecularRadiationMixedFvPatchScalarField:: +calcReceivedRayIDs() const +{ + if (receivedRayIDPtr_) + { + FatalErrorIn + ( + "wideBandSpecularRadiationMixedFvPatchScalarField" + "::calcreceivedRayIDs()" + ) << "receivedRayIDPtr already calculated" + << abort(FatalError); + } + + receivedRayIDPtr_ = new labelListList(this->size()); + labelListList& receivedRayID = *receivedRayIDPtr_; + + // Get access to radiation model, and recast as fvDOM + const radiationModel& radiation = + this->db().lookupObject("radiationProperties"); + + const fvDOM& dom(refCast(radiation)); + + // Get rayId and lambda Id for this ray + label rayId = -1; + label lambdaId = -1; + dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId); + + // Get index of active ray + word rayAndBand = this->dimensionedInternalField().name(); + rayAndBand = + rayAndBand.substr(rayAndBand.find("_")+1, rayAndBand.size()-1); + + // Get all ray directions + List dAve(dom.nRay()); + forAll (dAve, rayI) + { + dAve[rayI] = dom.IRay(rayI).dAve()/mag(dom.IRay(rayI).dAve()); + } + + // Get face normal vectors + vectorField nHat = this->patch().nf(); + + + // For each face, and for each reflected ray, try to find another + // ray that is better suited to pick it up. + forAll(receivedRayID, faceI) + { + // Access the list of received rays for this face + labelList& receivedRayIDs = receivedRayID[faceI]; + + // Check whether ray is going into surface + // -> no reflection + if ((dAve[rayId] & nHat[faceI]) > 0) + { + receivedRayIDs.setSize(0); + continue; + } + + + // For each face, initialize list of picked up + // rays to maximum possible size + receivedRayIDs.setSize(dom.nRay() - 1); + + // Count actual number of receiving rays + label nReceiving = 0; + + forAll(dAve, incomingRayI) + { + // Calculate reflected ray direction for rayI + vector dReflected = + dAve[incomingRayI] + - 2*(dAve[incomingRayI] & nHat[faceI])*nHat[faceI]; + + // Get dot product with this ray + scalar dotProductThisRay = dAve[rayId] & dReflected; + + // Assume this ray is closest to the reflected ray + bool closest = true; + + // And look through all other rays to find a better suited one + forAll(dAve, receivingRayI) + { + if (receivingRayI == rayId) + { + // Do not compare this ray with itself + continue; + } + + scalar dotProductOtherRay = dAve[receivingRayI] & dReflected; + if (dotProductThisRay < dotProductOtherRay) + { + // If another ray is closer, stop search + closest = false; + break; + } + } + + if (closest) + { + // Could not find better suited ray, so this ray needs to + // pick it up. Add incoming ray to list for this face + receivedRayIDs[nReceiving] = incomingRayI; + nReceiving ++; + } + + } + + // Resize list of picked up rays for this face + receivedRayIDs.setSize(nReceiving); + } + + // Sanity check on last ray + if (rayId == (dom.nRay() - 1)) + { + + forAll(nHat, faceI) + { + label incomingRays = 0; + label pickedUpRays = 0; + + for (label rayI = 0; rayI < dom.nRay(); rayI++) + { + // Is this ray going into the wall? + if ((dAve[rayI] & nHat[faceI]) > 0) + { + incomingRays++; + } + + // How many rays are picked up by this ray? + const wideBandSpecularRadiationMixedFvPatchScalarField& rayBC = + refCast + < + const wideBandSpecularRadiationMixedFvPatchScalarField + > + ( + dom.IRayLambda + ( + rayI, + lambdaId + ).boundaryField()[patch().index()] + ); + + const labelListList& receivedRaysList = rayBC.receivedRayIDs(); + + pickedUpRays += receivedRaysList[faceI].size(); + } + + if (incomingRays != pickedUpRays) + { + FatalErrorIn + ( + "wideBandSpecularRadiationMixedFvPatchScalarField" + "::calcreceivedRayIDs()" + ) << "Sanity checked failed on patch " << patch().name() + << " face " << faceI << nl + << "number of rays with direction into wall: " + << incomingRays + << ", number of reflected rays picked up by outgoing rays: " + << pickedUpRays << nl + << abort(FatalError); + } + } + } + +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +wideBandSpecularRadiationMixedFvPatchScalarField:: +wideBandSpecularRadiationMixedFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + mixedFvPatchField(p, iF), + receivedRayIDPtr_(NULL) +{ + this->refValue() = 0; + this->refGrad() = 0; + this->valueFraction() = 0; +} + + +wideBandSpecularRadiationMixedFvPatchScalarField:: +wideBandSpecularRadiationMixedFvPatchScalarField +( + const wideBandSpecularRadiationMixedFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + mixedFvPatchField(ptf, p, iF, mapper), + receivedRayIDPtr_(NULL) +{} + + +wideBandSpecularRadiationMixedFvPatchScalarField:: +wideBandSpecularRadiationMixedFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + mixedFvPatchField(p, iF), + receivedRayIDPtr_(NULL) +{ + this->refValue() = 0; + this->refGrad() = 0; + this->valueFraction() = 0; + + if (dict.found("patchType")) + { + word& pType = const_cast(this->patchType()); + pType = word(dict.lookup("patchType")); +} + + if (dict.found("value")) + { + fvPatchField::operator= + ( + scalarField("value", dict, p.size()) + ); + } + else + { + fvPatchScalarField::operator=(this->refValue()); + } +} + + +wideBandSpecularRadiationMixedFvPatchScalarField:: +wideBandSpecularRadiationMixedFvPatchScalarField +( + const wideBandSpecularRadiationMixedFvPatchScalarField& ptf, + const DimensionedField& iF +) +: + mixedFvPatchField(ptf, iF), + receivedRayIDPtr_(NULL) +{} + + +wideBandSpecularRadiationMixedFvPatchScalarField:: +wideBandSpecularRadiationMixedFvPatchScalarField +( + const wideBandSpecularRadiationMixedFvPatchScalarField& ptf +) +: + mixedFvPatchField(ptf), + receivedRayIDPtr_(NULL) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +const Foam::labelListList& +wideBandSpecularRadiationMixedFvPatchScalarField::receivedRayIDs() const +{ + if (!receivedRayIDPtr_) + { + calcReceivedRayIDs(); + } + + return *receivedRayIDPtr_; +} + + +// Evaluate the field on the patch +void wideBandSpecularRadiationMixedFvPatchScalarField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + vectorField nHat = this->patch().nf(); + + // Get access to radiation model, and recast as fvDOM + const radiationModel& radiation = + this->db().lookupObject("radiationProperties"); + + const fvDOM& dom(refCast(radiation)); + + // Get rayId and lambda Id for this ray + label rayId = -1; + label lambdaId = -1; + dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId); + + const labelListList& receivedRayIDs = this->receivedRayIDs(); + + scalarField IValue(this->size(), 0); + + labelList faceCellIDs = this->patch().faceCells(); + + // Loop over all faces and add values from all received rays + forAll(receivedRayIDs, faceI) + { + if (receivedRayIDs[faceI].size() == 0) + { + // Ray goes into face -> act as zeroGradient + this->valueFraction()[faceI] = 0; + } + else + { + // Ray goes out of face -> act as fixedValue + this->valueFraction()[faceI] = 1; + + // Pick up all reflected rays + forAll(receivedRayIDs[faceI], receivedRayI) + { + // Get ray from object registry + IValue[faceI] += + dom.IRayLambda + ( + receivedRayIDs[faceI][receivedRayI], + lambdaId + ).internalField()[faceCellIDs[faceI]]; + } + } + } + + // Set value for patch + this->refValue() = IValue; + + scalarField::operator= + ( + this->valueFraction()*this->refValue() + + + (1.0 - this->valueFraction())* + ( + this->patchInternalField() + + this->refGrad()/this->patch().deltaCoeffs() + ) + ); + + mixedFvPatchField::updateCoeffs(); +} + + +void wideBandSpecularRadiationMixedFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + this->writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +void wideBandSpecularRadiationMixedFvPatchScalarField::operator= +( + const fvPatchScalarField& ptf +) +{ + fvPatchScalarField::operator= + ( + this->valueFraction()*this->refValue() + + (1 - this->valueFraction())*ptf + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField +( + fvPatchScalarField, + wideBandSpecularRadiationMixedFvPatchScalarField +); + + +} // End namespace radiation +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.H new file mode 100644 index 000000000..2dfd05cd8 --- /dev/null +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.H @@ -0,0 +1,162 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Class + Foam::wideBandSpecularRadiationMixedFvPatchScalarField + +Description + Foam::wideBandSpecularRadiationMixedFvPatchScalarField + +SourceFiles + wideBandSpecularRadiationMixedFvPatchScalarField.C + +Author + Dominik Christ, Wikki Ltd. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#ifndef wideBandSpecularRadiationMixedFvPatchScalarField_H +#define wideBandSpecularRadiationMixedFvPatchScalarField_H + +#include "mixedFvPatchField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace radiation +{ + +/*---------------------------------------------------------------------------*\ + Class wideBandSpecularRadiationMixedFvPatchScalarField Declaration +\*---------------------------------------------------------------------------*/ + +class wideBandSpecularRadiationMixedFvPatchScalarField +: + public mixedFvPatchField +{ + + //- For each face, ray index of ray in direction of reflection + mutable labelListList *receivedRayIDPtr_; + + //- Calculate IDs of rays in direction of reflection + void calcReceivedRayIDs() const; + +public: + + //- Runtime type information + TypeName("wideBandSpecularRadiation"); + + + // Constructors + + //- Construct from patch and internal field + wideBandSpecularRadiationMixedFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + wideBandSpecularRadiationMixedFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given wideBandSpecularRadiationMixedFvPatchScalarField onto a new patch + wideBandSpecularRadiationMixedFvPatchScalarField + ( + const wideBandSpecularRadiationMixedFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + wideBandSpecularRadiationMixedFvPatchScalarField + ( + const wideBandSpecularRadiationMixedFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp > clone() const + { + return tmp > + ( + new wideBandSpecularRadiationMixedFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + wideBandSpecularRadiationMixedFvPatchScalarField + ( + const wideBandSpecularRadiationMixedFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp > clone + ( + const DimensionedField& iF + ) const + { + return tmp > + ( + new wideBandSpecularRadiationMixedFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Access + + //- Return ray IDs of reflected rays + const labelListList& receivedRayIDs() const; + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; + + + // Member operators + + virtual void operator=(const fvPatchScalarField& pvf); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace radiation +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +#endif + +// ************************************************************************* // diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C index 4b76fda39..25dda14ca 100644 --- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C +++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C @@ -202,7 +202,7 @@ Foam::radiation::blackBodyEmission::EbDeltaLambdaT const Vector2D& band ) const { - tmp Eb + tmp tEb ( new volScalarField ( @@ -217,28 +217,32 @@ Foam::radiation::blackBodyEmission::EbDeltaLambdaT radiation::sigmaSB*pow4(T) ) ); + volScalarField& Eb = tEb(); - - if (band == Vector2D::one) + if (magSqr(band - Vector2D::one) > SMALL) // Multiple bands? { - return Eb; - } - else - { - forAll(T, i) + forAll(Eb.internalField(), cellI) { - scalar T1 = fLambdaT(band[1]*T[i]); - scalar T2 = fLambdaT(band[0]*T[i]); - dimensionedScalar fLambdaDelta - ( - "fLambdaDelta", - dimless, - T1 - T2 - ); - Eb()[i] = Eb()[i]*fLambdaDelta.value(); + scalar percentileUpper = fLambdaT(band[1]*T.internalField()[cellI]); + scalar percentileLower = fLambdaT(band[0]*T.internalField()[cellI]); + Eb.internalField()[cellI] *= (percentileUpper - percentileLower); + } + + forAll(Eb.boundaryField(), patchI) + { + forAll(Eb.boundaryField()[patchI], faceI) + { + scalar percentileUpper = + fLambdaT(band[1]*T.boundaryField()[patchI][faceI]); + scalar percentileLower = + fLambdaT(band[0]*T.boundaryField()[patchI][faceI]); + Eb.boundaryField()[patchI][faceI] *= + (percentileUpper - percentileLower); + } } - return Eb; } + + return tEb; } diff --git a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C index fd58c8fc9..880383e38 100644 --- a/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C +++ b/src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C @@ -25,10 +25,12 @@ License #include "fvDOM.H" #include "addToRunTimeSelectionTable.H" - +#include "fvm.H" #include "mathematicalConstants.H" #include "radiationConstants.H" +using namespace Foam::mathematicalConstant; + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam @@ -57,66 +59,27 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T) IOobject ( "G", - mesh_.time().timeName(), + mesh().time().timeName(), T.db(), IOobject::NO_READ, IOobject::AUTO_WRITE ), - mesh_, + mesh(), dimensionedScalar("G", dimMass/pow3(dimTime), 0.0) ), - Qr_ - ( - IOobject - ( - "Qr", - mesh_.time().timeName(), - T.db(), - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - mesh_, - dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0) - ), a_ ( IOobject ( "a", - mesh_.time().timeName(), + mesh().time().timeName(), T.db(), IOobject::NO_READ, IOobject::AUTO_WRITE ), - mesh_, + mesh(), dimensionedScalar("a", dimless/dimLength, 0.0) ), - e_ - ( - IOobject - ( - "e", - mesh_.time().timeName(), - T.db(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh_, - dimensionedScalar("a", dimless/dimLength, 0.0) - ), - E_ - ( - IOobject - ( - "E", - mesh_.time().timeName(), - T.db(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh_, - dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0) - ), nTheta_(readLabel(coeffs_.lookup("nTheta"))), nPhi_(readLabel(coeffs_.lookup("nPhi"))), nRay_(0), @@ -124,15 +87,20 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T) aLambda_(nLambda_), blackBody_(nLambda_, T), IRay_(0), + Qem_(nLambda_), + Qin_(nLambda_), convergence_(coeffs_.lookupOrDefault("convergence", 0.0)), - maxIter_(coeffs_.lookupOrDefault