diff --git a/src/thermophysicalModels/radiation/Make/files b/src/thermophysicalModels/radiation/Make/files index 08c27f94bd..c1cd7aaf42 100644 --- a/src/thermophysicalModels/radiation/Make/files +++ b/src/thermophysicalModels/radiation/Make/files @@ -69,6 +69,7 @@ derivedFvPatchFields/MarshakRadiationFixedTemperature/MarshakRadiationFixedTempe derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C +derivedFvPatchFields/specularRadiation/specularRadiationMixedFvPatchScalarField.C /* fvOptions */ fvOptions/radiation/radiation.C diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/specularRadiation/specularRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/specularRadiation/specularRadiationMixedFvPatchScalarField.C new file mode 100644 index 0000000000..d748553358 --- /dev/null +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/specularRadiation/specularRadiationMixedFvPatchScalarField.C @@ -0,0 +1,468 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2023 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "radiationModel.H" +#include "fvDOM.H" +#include "specularRadiationMixedFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "wedgePolyPatch.H" +#include "symmetryPlanePolyPatch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace radiation +{ + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +scalar specularRadiationMixedFvPatchScalarField::azimuthAngle +( + const vector& d +) const +{ + return sign(d.y())*Foam::acos(d.x()/Foam::sqrt(sqr(d.x()) + sqr(d.y()))); +} + + +scalar specularRadiationMixedFvPatchScalarField::polarAngle +( + const vector& d +) const +{ + return Foam::acos(d.z()/mag(d)); +} + + +tmp specularRadiationMixedFvPatchScalarField::interpolateI +( + const fvDOM& dom, + const label closestRayi +) const +{ + // Calculate the reflected ray for this ray (KE:p. 2) + const vector dAve(normalised(dom.IRay(rayID_).dAve())); + const vector dSpe(normalised(dAve - 2*(dAve & n_)*n_)); + + + // Fetch the number of polar and azimuthal segments + const label nPolar = dom.nTheta(); + const label nAzimuth = 4*dom.nPhi(); + + + // Find the neighbouring ray indices of the reflected ray in the east + // and west + // Go to the east + const label polari = std::floor(closestRayi/nAzimuth) + 1; + label eastRayi = closestRayi + 1; + if (eastRayi == nAzimuth*polari) + { + eastRayi -= nAzimuth; + } + // Go to the west + label westRayi = closestRayi - 1; + if (westRayi < nAzimuth*(polari - 1)) + { + westRayi += nAzimuth; + } + + // Find the ray index closest to the reflected ray in the azimuthal + // direction + label azimuthRayi = -1; + bool east = false; + const vector dEast(normalised(dom.IRay(eastRayi).dAve())); + const vector dWest(normalised(dom.IRay(westRayi).dAve())); + if (mag(dSpe - dEast) < mag(dSpe - dWest)) + { + azimuthRayi = eastRayi; + east = true; + } + else + { + azimuthRayi = westRayi; + } + + + // Find the neighbouring ray indices of the reflected ray in the north + // and south + // Go to the north + label northRayi = closestRayi - nAzimuth; + if (northRayi < 0) + { + // The pole is inside the polar segment - skip + northRayi = -1; + } + // Go to the south + label southRayi = closestRayi + nAzimuth; + if (southRayi > nPolar*nAzimuth-1) + { + // The pole is inside the polar segment - skip + southRayi = -1; + } + + + // Find the ray index closest to the reflected ray in the polar direction + label polarRayi = -1; + if (northRayi != -1 && southRayi != -1) + { + const vector dNorth(normalised(dom.IRay(northRayi).dAve())); + const vector dSouth(normalised(dom.IRay(southRayi).dAve())); + + if (mag(dSpe - dNorth) < mag(dSpe - dSouth)) + { + polarRayi = northRayi; + } + else + { + polarRayi = southRayi; + } + } + else if (northRayi != -1) + { + polarRayi = northRayi; + } + else if (southRayi != -1) + { + polarRayi = southRayi; + } + + + // Find the ray index neighbouring the azimuthal and polar neighbour rays + label cornerRayi = -1; + if (polarRayi != -1) + { + if (east) + { + cornerRayi = polarRayi + 1; + + if (!(cornerRayi % nAzimuth)) + { + cornerRayi -= nAzimuth; + } + } + else + { + cornerRayi = polarRayi - 1; + if (!(polarRayi % nAzimuth)) + { + cornerRayi += nAzimuth; + } + } + } + + + // Interpolate the ray intensity of the reflected complementary ray from + // the neighbouring rays + const label patchi = this->patch().index(); + auto tIc = tmp::New(this->patch().size(), Zero); + auto& Ic = tIc.ref(); + + + if (polarRayi == -1) + { + // Linear interpolation only in the azimuth direction + const vector d1(normalised(dom.IRay(closestRayi).dAve())); + const vector d2(normalised(dom.IRay(azimuthRayi).dAve())); + + const scalar phic = azimuthAngle(dSpe); + const scalar phi1 = azimuthAngle(d1); + const scalar phi2 = azimuthAngle(d2); + + const auto& I1 = + dom.IRayLambda(closestRayi, lambdaID_).boundaryField()[patchi]; + const auto& I2 = + dom.IRayLambda(azimuthRayi, lambdaID_).boundaryField()[patchi]; + + Ic = lerp(I1, I2, (phic - phi1)/(phi2 - phi1)); + } + else + { + // Bilinear interpolation in the azimuth and polar directions + const vector d1(normalised(dom.IRay(closestRayi).dAve())); + const vector d2(normalised(dom.IRay(azimuthRayi).dAve())); + const vector d3(normalised(dom.IRay(polarRayi).dAve())); + const vector d4(normalised(dom.IRay(cornerRayi).dAve())); + + const scalar phic = azimuthAngle(dSpe); + const scalar phi1 = azimuthAngle(d1); + const scalar phi2 = azimuthAngle(d2); + const scalar phi3 = azimuthAngle(d3); + const scalar phi4 = azimuthAngle(d4); + + const scalar thetac = polarAngle(dSpe); + const scalar theta1 = polarAngle(d1); + const scalar theta3 = polarAngle(d3); + + const auto& I1 = + dom.IRayLambda(closestRayi, lambdaID_).boundaryField()[patchi]; + const auto& I2 = + dom.IRayLambda(azimuthRayi, lambdaID_).boundaryField()[patchi]; + const auto& I3 = + dom.IRayLambda(polarRayi, lambdaID_).boundaryField()[patchi]; + const auto& I4 = + dom.IRayLambda(cornerRayi, lambdaID_).boundaryField()[patchi]; + + const scalarField Ia(lerp(I1, I2, (phic - phi1)/(phi2 - phi1))); + const scalarField Ib(lerp(I3, I4, (phic - phi3)/(phi4 - phi3))); + + Ic = lerp(Ia, Ib, (thetac - theta1)/(theta3 - theta1)); + } + + return tIc; +} + + +label specularRadiationMixedFvPatchScalarField::calcComplementaryRayID +( + const fvDOM& dom +) const +{ + vectorList dAve(dom.nRay()); + forAll(dAve, rayi) + { + dAve[rayi] = normalised(dom.IRay(rayi).dAve()); + } + + + // Check if the ray goes out of the domain, ie. no reflection + if ((dAve[rayID_] & n_) > 0) + { + return -1; + } + + + // Calculate the reflected ray direction for rayID (KE:p. 2) + const vector dSpe + ( + normalised(dAve[rayID_] - 2*(dAve[rayID_] & n_)*n_) + ); + + + // Find the closest ray to the reflected ray + label complementaryRayID = -1; + scalar dotProductThisRay = -GREAT; + forAll(dAve, i) + { + const scalar dotProductOtherRay = dAve[i] & dSpe; + + if (dotProductThisRay < dotProductOtherRay) + { + complementaryRayID = i; + dotProductThisRay = dotProductOtherRay; + } + } + + return complementaryRayID; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +specularRadiationMixedFvPatchScalarField:: +specularRadiationMixedFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(p, iF), + n_(), + rayID_(-1), + lambdaID_(-1), + interpolate_(false) +{ + this->refValue() = Zero; + this->refGrad() = Zero; + this->valueFraction() = Zero; +} + + +specularRadiationMixedFvPatchScalarField:: +specularRadiationMixedFvPatchScalarField +( + const specularRadiationMixedFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + mixedFvPatchScalarField(ptf, p, iF, mapper), + n_(ptf.n_), + rayID_(ptf.rayID_), + lambdaID_(ptf.lambdaID_), + interpolate_(ptf.interpolate_) +{} + + +specularRadiationMixedFvPatchScalarField:: +specularRadiationMixedFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + mixedFvPatchScalarField(p, iF), + n_(), + rayID_(-1), + lambdaID_(-1), + interpolate_(dict.getOrDefault("interpolate", false)) +{ + this->refValue() = Zero; + this->refGrad() = Zero; + this->valueFraction() = Zero; + + if (!this->readValueEntry(dict)) + { + fvPatchScalarField::operator=(this->refValue()); + } + + + if (isA(p.patch())) + { + const auto& wp = refCast(p.patch()); + n_ = wp.n(); + } + else if (isA(p.patch())) + { + const auto& sp = refCast(p.patch()); + n_ = sp.n(); + } + else + { + FatalErrorInFunction + << " specularRadiation boundary condition is limited to " + << "wedge or symmetry-plane geometries." << nl + << exit(FatalError); + } +} + + +specularRadiationMixedFvPatchScalarField:: +specularRadiationMixedFvPatchScalarField +( + const specularRadiationMixedFvPatchScalarField& ptf, + const DimensionedField& iF +) +: + mixedFvPatchScalarField(ptf, iF), + n_(ptf.n_), + rayID_(ptf.rayID_), + lambdaID_(ptf.lambdaID_), + interpolate_(ptf.interpolate_) +{} + + +specularRadiationMixedFvPatchScalarField:: +specularRadiationMixedFvPatchScalarField +( + const specularRadiationMixedFvPatchScalarField& ptf +) +: + mixedFvPatchScalarField(ptf), + n_(ptf.n_), + rayID_(ptf.rayID_), + lambdaID_(ptf.lambdaID_), + interpolate_(ptf.interpolate_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void specularRadiationMixedFvPatchScalarField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + const fvDOM& dom = db().lookupObject("radiationProperties"); + + // Get rayID and lambdaID for this ray + if (rayID_ == -1 && lambdaID_ == -1) + { + dom.setRayIdLambdaId(internalField().name(), rayID_, lambdaID_); + } + + + // Find the ray index closest to the reflected ray + const label compID = calcComplementaryRayID(dom); + + + if (compID == -1) + { + // Apply zero-gradient condition for rays outgoing from the domain + this->valueFraction() = 0; + } + else + { + // Apply fixed condition for rays incoming to the domain + this->valueFraction() = 1; + + if (!interpolate_) + { + // Fetch the existing ray closest to the reflected ray (KE:Eq. 4) + this->refValue() = + dom.IRayLambda + ( + compID, + lambdaID_ + ).internalField(); + } + else + { + // Interpolate the ray intensity from neighbouring rays (KE:p. 2) + this->refValue() = interpolateI(dom, compID); + } + } + + mixedFvPatchScalarField::updateCoeffs(); +} + + +void specularRadiationMixedFvPatchScalarField::write(Ostream& os) const +{ + mixedFvPatchScalarField::write(os); + os.writeEntryIfDifferent("interpolate", false, interpolate_); + this->writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField +( + fvPatchScalarField, + specularRadiationMixedFvPatchScalarField +); + + +} // End namespace radiation +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/specularRadiation/specularRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiation/derivedFvPatchFields/specularRadiation/specularRadiationMixedFvPatchScalarField.H new file mode 100644 index 0000000000..bc456ae4fe --- /dev/null +++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/specularRadiation/specularRadiationMixedFvPatchScalarField.H @@ -0,0 +1,221 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2023 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::radiation::specularRadiationMixedFvPatchScalarField + +Description + This boundary condition provides a specular radiation condition for + axisymmetric and symmetry-plane \c fvDOM computations. + + References: + \verbatim + Standard model (tag:KE): + Kumar, P., & Eswaran, V. (2013). + A methodology to solve 2D and axisymmetric radiative + transfer problems using a general 3D solver. + Journal of heat transfer, 135(12). + DOI:10.1115/1.4024674 + \endverbatim + +Usage + Example of the boundary condition specification: + \verbatim + + { + // Mandatory entries + type specularRadiation; + + // Optional entries + interpolate ; + + // Inherited entries + patchType ; + ... + } + \endverbatim + + where the entries mean: + \table + Property | Description | Type | Reqd | Deflt + type | Type name: turbulentDigitalFilterInlet | word | yes | - + interpolate | Flag to enable ray-intensity interp | bool | no | false + \endtable + + The inherited entries are elaborated in: + - \link mixedFvPatchFields.H \endlink + +Note + - The condition is limited to the underlying \c wedge and \c symmetryPlane + conditions. + +SourceFiles + specularRadiationMixedFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_specularRadiationMixedFvPatchScalarField_H +#define Foam_specularRadiationMixedFvPatchScalarField_H + +#include "mixedFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace radiation +{ + +/*---------------------------------------------------------------------------*\ + Class specularRadiationMixedFvPatchScalarField Declaration +\*---------------------------------------------------------------------------*/ + +class specularRadiationMixedFvPatchScalarField +: + public mixedFvPatchScalarField +{ + // Private Data + + //- Patch normal vector + vector n_; + + //- Ray index of this ray + label rayID_; + + //- Band index of this ray + label lambdaID_; + + //- Flag to enable ray-intensity interpolation + const bool interpolate_; + + + // Private Member Functions + + //- Return the corresponding azimuth angle of a given Cartesian vector + scalar azimuthAngle(const vector& d) const; + + //- Return the corresponding polar angle of a given Cartesian vector + scalar polarAngle(const vector& d) const; + + //- Return interpolated ray intensity + tmp interpolateI + ( + const fvDOM& dom, + const label closestRayID + ) const; + + //- Return the index of complementary ray + label calcComplementaryRayID(const fvDOM& dom) const; + + +public: + + //- Runtime type information + TypeName("specularRadiation"); + + + // Constructors + + //- Construct from patch and internal field + specularRadiationMixedFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + specularRadiationMixedFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given + //- specularRadiationMixedFvPatchScalarField onto a new patch + specularRadiationMixedFvPatchScalarField + ( + const specularRadiationMixedFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + specularRadiationMixedFvPatchScalarField + ( + const specularRadiationMixedFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp > clone() const + { + return tmp > + ( + new specularRadiationMixedFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + specularRadiationMixedFvPatchScalarField + ( + const specularRadiationMixedFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp > clone + ( + const DimensionedField& iF + ) const + { + return tmp > + ( + new specularRadiationMixedFvPatchScalarField(*this, iF) + ); + } + + + // Member Functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace radiation +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +#endif + +// ************************************************************************* //