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;
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
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;
615 List<bool> rayIdConv(nRay_,
false);
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_();
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)
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.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
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)
Internal & ref(const bool updateAccessTime=true)
Same as internalFieldRef()
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)