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
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;
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;
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_(),
450 Foam::radiation::fvDOM::fvDOM
462 mesh_.time().timeName(),
475 mesh_.time().timeName(),
488 mesh_.time().timeName(),
501 mesh_.time().timeName(),
514 mesh_.time().timeName(),
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)()();
733 void 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_();