ENH: specularRadiation: add axisymmetric fvDOM boundary condition

This commit is contained in:
Kutalmis Bercin
2023-06-02 13:54:02 +01:00
committed by Andrew Heather
parent 8f46e47931
commit ee26a36add
3 changed files with 690 additions and 0 deletions

View File

@ -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

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<scalarField> 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<scalarField>::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<scalar, volMesh>& 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<scalar, volMesh>& 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<scalar, volMesh>& 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<wedgePolyPatch>(p.patch()))
{
const auto& wp = refCast<const wedgePolyPatch>(p.patch());
n_ = wp.n();
}
else if (isA<symmetryPlanePolyPatch>(p.patch()))
{
const auto& sp = refCast<const symmetryPlanePolyPatch>(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<scalar, volMesh>& 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<fvDOM>("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<bool>("interpolate", false, interpolate_);
this->writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchScalarField,
specularRadiationMixedFvPatchScalarField
);
} // End namespace radiation
} // End namespace Foam
// ************************************************************************* //

View File

@ -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 <http://www.gnu.org/licenses/>.
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
<patchName>
{
// Mandatory entries
type specularRadiation;
// Optional entries
interpolate <bool>;
// Inherited entries
patchType <word>;
...
}
\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<scalarField> 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<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
specularRadiationMixedFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given
//- specularRadiationMixedFvPatchScalarField onto a new patch
specularRadiationMixedFvPatchScalarField
(
const specularRadiationMixedFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
specularRadiationMixedFvPatchScalarField
(
const specularRadiationMixedFvPatchScalarField&
);
//- Construct and return a clone
virtual tmp<fvPatchField<scalar> > clone() const
{
return tmp<fvPatchField<scalar> >
(
new specularRadiationMixedFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
specularRadiationMixedFvPatchScalarField
(
const specularRadiationMixedFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchField<scalar> > clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchField<scalar> >
(
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
// ************************************************************************* //