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);
117void Foam::radiation::fvDOM::initialise()
119 coeffs_.readIfPresent(
"useExternalBeam", useExternalBeam_);
121 if (useExternalBeam_)
123 spectralDistributions_.reset
127 "spectralDistribution",
133 spectralDistribution_ =
134 spectralDistributions_->value(mesh_.time().timeOutputValue());
136 spectralDistribution_ =
137 spectralDistribution_/
sum(spectralDistribution_);
139 const dictionary& solarDict = this->subDict(
"solarCalculatorCoeffs");
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;
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;
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;
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;
416 nTheta_(coeffs_.get<label>(
"nTheta")),
417 nPhi_(coeffs_.get<label>(
"nPhi")),
419 nLambda_(absorptionEmission_->nBands()),
421 blackBody_(nLambda_,
T),
425 coeffs_.getOrDefaultCompat<scalar>
428 {{
"convergence", 1712}},
432 maxIter_(coeffs_.getOrDefault<label>(
"maxIter", 50)),
434 useSolarLoad_(
false),
438 coeffs_.getOrDefault<
vector>(
"meshOrientation",
Zero)
440 useExternalBeam_(
false),
441 spectralDistribution_(),
442 spectralDistributions_(),
522 nTheta_(coeffs_.get<label>(
"nTheta")),
523 nPhi_(coeffs_.get<label>(
"nPhi")),
525 nLambda_(absorptionEmission_->nBands()),
527 blackBody_(nLambda_,
T),
531 coeffs_.getOrDefaultCompat<scalar>
534 {{
"convergence", 1712}},
538 maxIter_(coeffs_.getOrDefault<label>(
"maxIter", 50)),
540 useSolarLoad_(
false),
544 coeffs_.getOrDefault<
vector>(
"meshOrientation",
Zero)
546 useExternalBeam_(
false),
547 spectralDistribution_(),
548 spectralDistributions_(),
563 coeffs_.readIfPresentCompat
565 "tolerance", {{
"convergence", 1712}}, tolerance_
567 coeffs_.readIfPresent(
"maxIter", maxIter_);
578 absorptionEmission_->correct(a_, aLambda_);
580 updateBlackBodyEmission();
584 solarLoad_->calculate();
587 if (useExternalBeam_)
589 switch (solarCalculator_->sunDirectionModel())
597 label updateIndex = label
600 /solarCalculator_->sunTrackingUpdateInterval()
603 if (updateIndex > updateTimeIndex_)
605 Info <<
"Updating Sun position..." <<
endl;
606 updateTimeIndex_ = updateIndex;
617 scalar maxResidual = 0;
621 Info<<
"Radiation solver iter: " << radIter <<
endl;
627 if (!rayIdConv[rayI])
629 scalar maxBandResidual = IRay_[rayI].correct();
630 maxResidual =
max(maxBandResidual, maxResidual);
632 if (maxBandResidual < tolerance_)
634 rayIdConv[rayI] =
true;
639 }
while (maxResidual > tolerance_ && radIter < maxIter_);
655 mesh_.time().timeName(),
664 *(aLambda_[0] - absorptionEmission_->aDisp(0)())
665 *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(0))
673 for (label j=1; j < nLambda_; j++)
679 *(aLambda_[j] - absorptionEmission_->aDisp(j)())
680 *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(j))
698 mesh_.time().timeName(),
712 for (label j=0; j < nLambda_; j++)
717 IRay_[0].ILambda(j)()*IRay_[0].omega()
720 for (label rayI=1; rayI < nRay_; rayI++)
722 Gj.
ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
725 Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
726 - absorptionEmission_->ECont(j)()();
733void Foam::radiation::fvDOM::updateBlackBodyEmission()
735 for (label j=0; j < nLambda_; j++)
737 blackBody_.correct(j, absorptionEmission_->bands(j));
751 IRay_[rayI].addIntensity();
752 G_ += IRay_[rayI].I()*IRay_[rayI].omega();
753 qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
754 qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
755 qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
768 const auto i1 =
name.find(
'_');
769 const auto i2 =
name.find(
'_', i1+1);
778 return solarCalculator_();
Macros for easy insertion into run-time selection tables.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
Finite Volume Discrete Ordinates Method. Solves the RTE equation for n directions in a participating ...
virtual tmp< volScalarField::Internal > Ru() const
Source term component (constant)
virtual tmp< volScalarField > Rp() const
Source term component (for power of T^4)
const solarCalculator & solarCalc() const
Solar calculator.
void alignClosestRayToSun(const vector &sunDir)
Align closest ray to sunDir.
void updateG()
Update G and calculate total heat flux on boundary.
void rotateInitialRays(const vector &sunDir)
Rotate rays spheric equator to sunDir.
void updateRaysDir()
Rotate rays according to Sun direction.
void setRayIdLambdaId(const word &name, label &rayId, label &lambdaId) const
Set the rayId and lambdaId from by decomposing an intensity.
bool read()
Read radiation properties dictionary.
void calculate()
Solve radiation equation(s)
Top level model for radiation modelling.
virtual bool read()=0
Read radiationProperties dictionary.
Radiation intensity for a ray in a given direction.
The solarLoad radiation model includes Sun primary hits, their reflective fluxes and diffusive sky ra...
A solar calculator model providing models for the solar direction and solar loads.
A class for managing temporary objects.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
autoPtr< radiation::radiationModel > radiation(radiation::radiationModel::New(T))
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
constexpr scalar pi(M_PI)
constexpr scalar piByTwo(0.5 *M_PI)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
Different types of constants.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
const dimensionSet dimless
Dimensionless.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensionedScalar pow3(const dimensionedScalar &ds)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
label readLabel(const char *buf)
Parse entire buffer as a label, skipping leading/trailing whitespace.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
errorManip< error > abort(error &err)
static constexpr const zero Zero
Global zero (0)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
constexpr char nl
The newline '\n' character (0x0a)
#define addToRadiationRunTimeSelectionTables(model)
#define forAll(list, i)
Loop across all elements in list.
Unit conversion functions.