62 IRay_[rayId].dAve() = coordRot & IRay_[rayId].dAve();
63 IRay_[rayId].d() = coordRot & IRay_[rayId].d();
71 scalar maxSunRay = -GREAT;
76 const vector& iD = IRay_[rayId].d();
77 scalar dir = sunDir & iD;
90 IRay_[rayId].dAve() = coordRot & IRay_[rayId].dAve();
91 IRay_[rayId].d() = coordRot & IRay_[rayId].d();
94 Info <<
"Sun direction : " << sunDir <<
nl <<
endl;
95 Info <<
"Sun ray ID : " << SunRayId <<
nl <<
endl;
101 solarCalculator_->correctSunDirection();
102 const vector sunDir = solarCalculator_->direction();
105 if (updateTimeIndex_ == 0)
107 rotateInitialRays(sunDir);
108 alignClosestRayToSun(sunDir);
110 else if (updateTimeIndex_ > 0)
112 alignClosestRayToSun(sunDir);
117 void Foam::radiation::fvDOM::initialise()
119 coeffs_.readIfPresent(
"useExternalBeam", useExternalBeam_);
121 if (useExternalBeam_)
123 spectralDistributions_.reset
125 Function1<scalarField>::New
127 "spectralDistribution",
133 spectralDistribution_ =
134 spectralDistributions_->value(mesh_.time().timeOutputValue());
136 spectralDistribution_ =
137 spectralDistribution_/
sum(spectralDistribution_);
139 const dictionary& solarDict = this->subDict(
"solarCalculatorCoeffs");
140 solarCalculator_.reset(
new solarCalculator(solarDict, mesh_));
142 if (mesh_.nSolutionD() != 3)
145 <<
"External beam model only available in 3D meshes " 149 if (solarCalculator_->diffuseSolarRad() > 0)
152 <<
"External beam model does not support Diffuse " 153 <<
"Solar Radiation. Set diffuseSolarRad to zero" 156 if (spectralDistribution_.size() != nLambda_)
159 <<
"The epectral energy distribution has different bands " 160 <<
"than the absoprtivity model " 166 if (mesh_.nSolutionD() == 3)
168 nRay_ = 4*nPhi_*nTheta_;
170 IRay_.setSize(nRay_);
172 const scalar deltaPhi =
pi/(2*nPhi_);
173 const scalar deltaTheta =
pi/nTheta_;
177 for (label
n = 1;
n <= nTheta_;
n++)
179 for (label m = 1; m <= 4*nPhi_; m++)
181 scalar thetai = (2*
n - 1)*deltaTheta/2.0;
182 scalar phii = (2*m - 1)*deltaPhi/2.0;
187 new radiativeIntensityRay
196 *absorptionEmission_,
206 else if (mesh_.nSolutionD() == 2)
209 const scalar deltaTheta =
pi;
211 IRay_.setSize(nRay_);
212 const scalar deltaPhi =
pi/(2.0*nPhi_);
214 for (label m = 1; m <= 4*nPhi_; m++)
216 const scalar phii = (2*m - 1)*deltaPhi/2.0;
220 new radiativeIntensityRay
229 *absorptionEmission_,
241 const scalar deltaTheta =
pi;
243 IRay_.setSize(nRay_);
244 const scalar deltaPhi =
pi;
246 for (label m = 1; m <= 2; m++)
248 const scalar phii = (2*m - 1)*deltaPhi/2.0;
252 new radiativeIntensityRay
261 *absorptionEmission_,
282 mesh_.time().timeName(),
292 Info<<
"fvDOM : Allocated " << IRay_.size()
293 <<
" rays with average orientation:" <<
nl;
295 if (useExternalBeam_)
301 scalar totalOmega = 0;
304 if (omegaMax_ < IRay_[rayId].omega())
306 omegaMax_ = IRay_[rayId].omega();
308 totalOmega += IRay_[rayId].omega();
309 Info<<
'\t' << IRay_[rayId].I().name() <<
" : " <<
"dAve : " 310 <<
'\t' << IRay_[rayId].dAve() <<
" : " <<
"omega : " 311 <<
'\t' << IRay_[rayId].omega() <<
" : " <<
"d : " 312 <<
'\t' << IRay_[rayId].d() <<
nl;
315 Info <<
"Total omega : " << totalOmega <<
endl;
319 coeffs_.readIfPresent(
"useSolarLoad", useSolarLoad_);
323 if (useExternalBeam_)
326 <<
"External beam with fvDOM can not be used " 327 <<
"with the solar load model" 330 const dictionary& solarDict = this->subDict(
"solarLoadCoeffs");
331 solarLoad_.reset(
new solarLoad(solarDict, T_));
333 if (solarLoad_->nBands() != this->nBands())
336 <<
"Requested solar radiation with fvDOM. Using " 337 <<
"different number of bands for the solar load is not allowed" 341 Info<<
"Creating Solar Load Model " <<
nl;
419 nTheta_(coeffs_.
get<label>(
"nTheta")),
420 nPhi_(coeffs_.
get<label>(
"nPhi")),
422 nLambda_(absorptionEmission_->nBands()),
424 blackBody_(nLambda_,
T),
428 coeffs_.getOrDefaultCompat<scalar>
431 {{
"convergence", 1712}},
435 maxIter_(coeffs_.getOrDefault<label>(
"maxIter", 50)),
437 useSolarLoad_(
false),
441 coeffs_.getOrDefault<
vector>(
"meshOrientation",
Zero)
443 useExternalBeam_(
false),
444 spectralDistribution_(),
445 spectralDistributions_(),
453 Foam::radiation::fvDOM::fvDOM
528 nTheta_(coeffs_.
get<label>(
"nTheta")),
529 nPhi_(coeffs_.
get<label>(
"nPhi")),
531 nLambda_(absorptionEmission_->nBands()),
533 blackBody_(nLambda_,
T),
537 coeffs_.getOrDefaultCompat<scalar>
540 {{
"convergence", 1712}},
544 maxIter_(coeffs_.getOrDefault<label>(
"maxIter", 50)),
546 useSolarLoad_(
false),
550 coeffs_.getOrDefault<
vector>(
"meshOrientation",
Zero)
552 useExternalBeam_(
false),
553 spectralDistribution_(),
554 spectralDistributions_(),
569 coeffs_.readIfPresentCompat
571 "tolerance", {{
"convergence", 1712}}, tolerance_
573 coeffs_.readIfPresent(
"maxIter", maxIter_);
584 absorptionEmission_->correct(a_, aLambda_);
586 updateBlackBodyEmission();
590 solarLoad_->calculate();
593 if (useExternalBeam_)
595 switch (solarCalculator_->sunDirectionModel())
603 label updateIndex = label
606 /solarCalculator_->sunTrackingUpdateInterval()
609 if (updateIndex > updateTimeIndex_)
611 Info <<
"Updating Sun position..." <<
endl;
612 updateTimeIndex_ = updateIndex;
621 List<bool> rayIdConv(nRay_,
false);
623 scalar maxResidual = 0;
627 Info<<
"Radiation solver iter: " << radIter <<
endl;
633 if (!rayIdConv[rayI])
635 scalar maxBandResidual = IRay_[rayI].correct();
636 maxResidual =
max(maxBandResidual, maxResidual);
638 if (maxBandResidual < tolerance_)
640 rayIdConv[rayI] =
true;
645 }
while (maxResidual > tolerance_ && radIter < maxIter_);
661 * (aLambda_[0] - absorptionEmission_->aDisp(0)())
662 * blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(0))
665 auto& Rp = tRp.ref();
668 for (label j=1; j < nLambda_; j++)
674 *(aLambda_[j] - absorptionEmission_->aDisp(j)())
675 *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(j))
693 auto& Ru = tRu.ref();
696 for (label j=0; j < nLambda_; j++)
701 IRay_[0].ILambda(j)()*IRay_[0].omega()
704 for (label rayI=1; rayI < nRay_; rayI++)
706 Gj.
ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
709 Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
710 - absorptionEmission_->ECont(j)()();
717 void Foam::radiation::fvDOM::updateBlackBodyEmission()
719 for (label j=0; j < nLambda_; j++)
721 blackBody_.correct(j, absorptionEmission_->bands(j));
735 IRay_[rayI].addIntensity();
736 G_ += IRay_[rayI].I()*IRay_[rayI].omega();
737 qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
738 qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
739 qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
752 const auto i1 =
name.find(
'_');
753 const auto i2 =
name.find(
'_', i1+1);
762 return solarCalculator_();
Different types of constants.
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
A solar calculator model providing models for the solar direction and solar loads.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Unit conversion functions.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
label readLabel(const char *buf)
Parse entire buffer as a label, skipping leading/trailing whitespace.
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
void setRayIdLambdaId(const word &name, label &rayId, label &lambdaId) const
Set the rayId and lambdaId from by decomposing an intensity.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
virtual bool read()=0
Read radiationProperties dictionary.
A class for handling words, derived from Foam::string.
void updateG()
Update G and calculate total heat flux on boundary.
constexpr scalar pi(M_PI)
Top level model for radiation modelling.
void rotateInitialRays(const vector &sunDir)
Rotate rays spheric equator to sunDir.
constexpr scalar piByTwo(0.5 *M_PI)
errorManip< error > abort(error &err)
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions and patch type. [Takes current timeN...
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
bool read()
Read radiation properties dictionary.
defineTypeNameAndDebug(combustionModel, 0)
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
void updateRaysDir()
Rotate rays according to Sun direction.
dimensionedScalar pow3(const dimensionedScalar &ds)
virtual tmp< volScalarField::Internal > Ru() const
Source term component (constant)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const Field< Type > &iField)
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions, copy of internal field...
const solarCalculator & solarCalc() const
Solar calculator.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
void calculate()
Solve radiation equation(s)
messageStream Info
Information stream (stdout output on master, null elsewhere)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
A class for managing temporary objects.
Tensor of scalars, i.e. Tensor<scalar>.
Defines the attributes of an object for which implicit objectRegistry management is supported...
virtual tmp< volScalarField > Rp() const
Source term component (for power of T^4)
#define addToRadiationRunTimeSelectionTables(model)
Finite Volume Discrete Ordinates Method. Solves the RTE equation for n directions in a participating ...
Do not request registration (bool: false)
void alignClosestRayToSun(const vector &sunDir)
Align closest ray to sunDir.
autoPtr< radiation::radiationModel > radiation(radiation::radiationModel::New(T))
static constexpr const zero Zero
Global zero (0)