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 coeffs_.readEntry(
"spectralDistribution", spectralDistribution_);
125 spectralDistribution_ =
126 spectralDistribution_/
sum(spectralDistribution_);
128 const dictionary& solarDict = this->subDict(
"solarCalculatorCoeffs");
131 if (mesh_.nSolutionD() != 3)
134 <<
"External beam model only available in 3D meshes "
138 if (solarCalculator_->diffuseSolarRad() > 0)
141 <<
"External beam model does not support Diffuse "
142 <<
"Solar Radiation. Set diffuseSolarRad to zero"
145 if (spectralDistribution_.size() != nLambda_)
148 <<
"The epectral energy distribution has different bands "
149 <<
"than the absoprtivity model "
155 if (mesh_.nSolutionD() == 3)
157 nRay_ = 4*nPhi_*nTheta_;
159 IRay_.setSize(nRay_);
161 const scalar deltaPhi =
pi/(2*nPhi_);
162 const scalar deltaTheta =
pi/nTheta_;
166 for (label
n = 1;
n <= nTheta_;
n++)
168 for (label m = 1; m <= 4*nPhi_; m++)
170 scalar thetai = (2*
n - 1)*deltaTheta/2.0;
171 scalar phii = (2*m - 1)*deltaPhi/2.0;
176 new radiativeIntensityRay
185 *absorptionEmission_,
195 else if (mesh_.nSolutionD() == 2)
198 const scalar deltaTheta =
pi;
200 IRay_.setSize(nRay_);
201 const scalar deltaPhi =
pi/(2.0*nPhi_);
203 for (label m = 1; m <= 4*nPhi_; m++)
205 const scalar phii = (2*m - 1)*deltaPhi/2.0;
209 new radiativeIntensityRay
218 *absorptionEmission_,
230 const scalar deltaTheta =
pi;
232 IRay_.setSize(nRay_);
233 const scalar deltaPhi =
pi;
235 for (label m = 1; m <= 2; m++)
237 const scalar phii = (2*m - 1)*deltaPhi/2.0;
241 new radiativeIntensityRay
250 *absorptionEmission_,
271 mesh_.time().timeName(),
281 Info<<
"fvDOM : Allocated " << IRay_.size()
282 <<
" rays with average orientation:" <<
nl;
284 if (useExternalBeam_)
290 scalar totalOmega = 0;
293 if (omegaMax_ < IRay_[rayId].omega())
295 omegaMax_ = IRay_[rayId].omega();
297 totalOmega += IRay_[rayId].omega();
298 Info<<
'\t' << IRay_[rayId].I().name() <<
" : " <<
"dAve : "
299 <<
'\t' << IRay_[rayId].dAve() <<
" : " <<
"omega : "
300 <<
'\t' << IRay_[rayId].omega() <<
" : " <<
"d : "
301 <<
'\t' << IRay_[rayId].d() <<
nl;
304 Info <<
"Total omega : " << totalOmega <<
endl;
308 coeffs_.readIfPresent(
"useSolarLoad", useSolarLoad_);
312 if (useExternalBeam_)
315 <<
"External beam with fvDOM can not be used "
316 <<
"with the solar load model"
319 const dictionary& solarDict = this->subDict(
"solarLoadCoeffs");
320 solarLoad_.reset(
new solarLoad(solarDict, T_));
322 if (solarLoad_->nBands() != this->nBands())
325 <<
"Requested solar radiation with fvDOM. Using "
326 <<
"different number of bands for the solar load is not allowed"
330 Info<<
"Creating Solar Load Model " <<
nl;
405 nTheta_(coeffs_.get<label>(
"nTheta")),
406 nPhi_(coeffs_.get<label>(
"nPhi")),
408 nLambda_(absorptionEmission_->nBands()),
410 blackBody_(nLambda_,
T),
414 coeffs_.getOrDefaultCompat<scalar>
417 {{
"convergence", 1712}},
421 maxIter_(coeffs_.getOrDefault<label>(
"maxIter", 50)),
423 useSolarLoad_(
false),
427 coeffs_.getOrDefault<
vector>(
"meshOrientation",
Zero)
429 useExternalBeam_(
false),
430 spectralDistribution_(),
438 Foam::radiation::fvDOM::fvDOM
450 mesh_.time().timeName(),
463 mesh_.time().timeName(),
476 mesh_.time().timeName(),
489 mesh_.time().timeName(),
502 mesh_.time().timeName(),
510 nTheta_(coeffs_.get<label>(
"nTheta")),
511 nPhi_(coeffs_.get<label>(
"nPhi")),
513 nLambda_(absorptionEmission_->nBands()),
515 blackBody_(nLambda_,
T),
519 coeffs_.getOrDefaultCompat<scalar>
522 {{
"convergence", 1712}},
526 maxIter_(coeffs_.getOrDefault<label>(
"maxIter", 50)),
528 useSolarLoad_(
false),
532 coeffs_.getOrDefault<
vector>(
"meshOrientation",
Zero)
534 useExternalBeam_(
false),
535 spectralDistribution_(),
556 coeffs_.readIfPresentCompat
558 "tolerance", {{
"convergence", 1712}}, tolerance_
560 coeffs_.readIfPresent(
"maxIter", maxIter_);
571 absorptionEmission_->correct(a_, aLambda_);
573 updateBlackBodyEmission();
577 solarLoad_->calculate();
580 if (useExternalBeam_)
582 switch (solarCalculator_->sunDirectionModel())
590 label updateIndex = label
593 /solarCalculator_->sunTrackingUpdateInterval()
596 if (updateIndex > updateTimeIndex_)
598 Info <<
"Updating Sun position..." <<
endl;
599 updateTimeIndex_ = updateIndex;
610 scalar maxResidual = 0;
614 Info<<
"Radiation solver iter: " << radIter <<
endl;
620 if (!rayIdConv[rayI])
622 scalar maxBandResidual = IRay_[rayI].correct();
623 maxResidual =
max(maxBandResidual, maxResidual);
625 if (maxBandResidual < tolerance_)
627 rayIdConv[rayI] =
true;
632 }
while (maxResidual > tolerance_ && radIter < maxIter_);
648 mesh_.time().timeName(),
657 *(aLambda_[0] - absorptionEmission_->aDisp(0)())
658 *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(0))
666 for (label j=1; j < nLambda_; j++)
672 *(aLambda_[j] - absorptionEmission_->aDisp(j)())
673 *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(j))
691 mesh_.time().timeName(),
705 for (label j=0; j < nLambda_; j++)
710 IRay_[0].ILambda(j)()*IRay_[0].omega()
713 for (label rayI=1; rayI < nRay_; rayI++)
715 Gj.
ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
718 Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
719 - absorptionEmission_->ECont(j)()();
726 void Foam::radiation::fvDOM::updateBlackBodyEmission()
728 for (label j=0; j < nLambda_; j++)
730 blackBody_.correct(j, absorptionEmission_->bands(j));
744 IRay_[rayI].addIntensity();
745 G_ += IRay_[rayI].I()*IRay_[rayI].omega();
746 qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
747 qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
748 qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
761 const auto i1 =
name.find(
'_');
762 const auto i2 =
name.find(
'_', i1+1);
771 return solarCalculator_();