fvDOM.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "fvDOM.H"
28 #include "scatterModel.H"
29 #include "constants.H"
30 #include "fvm.H"
31 #include "wedgePolyPatch.H"
33 
34 using namespace Foam::constant;
35 using namespace Foam::constant::mathematical;
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 namespace radiationModels
42 {
43  defineTypeNameAndDebug(fvDOM, 0);
45 }
46 }
47 
48 
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 
51 void Foam::radiationModels::fvDOM::initialise()
52 {
53  forAll(mesh_.boundaryMesh(), patchi)
54  {
55  const polyPatch& pp = mesh_.boundaryMesh()[patchi];
56  if
57  (
58  (
59  isA<cyclicTransform>(pp)
60  && refCast<const cyclicTransform>(pp).transform().rotates()
61  )
62  || isA<wedgePolyPatch>(pp)
63  )
64  {
66  << " radiation model does not currently support"
67  " rotationally transforming patches: cyclic and wedge."
68  << exit(FatalError);
69  }
70  }
71 
72  // 3D
73  if (mesh_.nSolutionD() == 3)
74  {
75  nRay_ = 4*nPhi_*nTheta_;
76  IRay_.setSize(nRay_);
77  const scalar deltaPhi = pi/(2.0*nPhi_);
78  const scalar deltaTheta = pi/nTheta_;
79  label i = 0;
80  for (label n = 1; n <= nTheta_; n++)
81  {
82  for (label m = 1; m <= 4*nPhi_; m++)
83  {
84  const scalar thetai = (2*n - 1)*deltaTheta/2.0;
85  const scalar phii = (2*m - 1)*deltaPhi/2.0;
86  IRay_.set
87  (
88  i,
89  new radiativeIntensityRay
90  (
91  *this,
92  mesh_,
93  phii,
94  thetai,
95  deltaPhi,
96  deltaTheta,
97  nLambda_,
98  absorptionEmission_,
99  blackBody_,
100  i
101  )
102  );
103  i++;
104  }
105  }
106  }
107  // 2D
108  else if (mesh_.nSolutionD() == 2)
109  {
110  const scalar thetai = piByTwo;
111  const scalar deltaTheta = pi;
112  nRay_ = 4*nPhi_;
113  IRay_.setSize(nRay_);
114  const scalar deltaPhi = pi/(2.0*nPhi_);
115  label i = 0;
116  for (label m = 1; m <= 4*nPhi_; m++)
117  {
118  const scalar phii = (2*m - 1)*deltaPhi/2.0;
119  IRay_.set
120  (
121  i,
122  new radiativeIntensityRay
123  (
124  *this,
125  mesh_,
126  phii,
127  thetai,
128  deltaPhi,
129  deltaTheta,
130  nLambda_,
131  absorptionEmission_,
132  blackBody_,
133  i
134  )
135  );
136  i++;
137  }
138  }
139  // 1D
140  else
141  {
142  const scalar thetai = piByTwo;
143  const scalar deltaTheta = pi;
144  nRay_ = 2;
145  IRay_.setSize(nRay_);
146  const scalar deltaPhi = pi;
147  label i = 0;
148  for (label m = 1; m <= 2; m++)
149  {
150  const scalar phii = (2*m - 1)*deltaPhi/2.0;
151  IRay_.set
152  (
153  i,
154  new radiativeIntensityRay
155  (
156  *this,
157  mesh_,
158  phii,
159  thetai,
160  deltaPhi,
161  deltaTheta,
162  nLambda_,
163  absorptionEmission_,
164  blackBody_,
165  i
166  )
167  );
168  i++;
169  }
170  }
171 
172 
173  // Construct absorption field for each wavelength
174  forAll(aLambda_, lambdaI)
175  {
176  aLambda_.set
177  (
178  lambdaI,
179  new volScalarField
180  (
181  IOobject
182  (
183  "aLambda_" + Foam::name(lambdaI) ,
184  mesh_.time().timeName(),
185  mesh_,
188  ),
189  a_
190  )
191  );
192  }
193 
194 
195  // Calculate the maximum solid angle
196  forAll(IRay_, rayId)
197  {
198  if (omegaMax_ < IRay_[rayId].omega())
199  {
200  omegaMax_ = IRay_[rayId].omega();
201  }
202  }
203 
204 
205  Info<< typeName << ": Created " << IRay_.size() << " rays with average "
206  << "directions (dAve) and solid angles (omega)" << endl;
207  Info<< incrIndent;
208  forAll(IRay_, rayId)
209  {
210  Info<< indent
211  << "Ray " << IRay_[rayId].I().name() << ": "
212  << "dAve = " << IRay_[rayId].dAve() << ", "
213  << "omega = " << IRay_[rayId].omega() << endl;
214  }
215  Info<< decrIndent << endl;
216 }
217 
218 
219 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
220 
222 :
223  radiationModel(typeName, T),
224  G_
225  (
226  IOobject
227  (
228  "G",
229  mesh_.time().timeName(),
230  mesh_,
231  IOobject::NO_READ,
232  IOobject::AUTO_WRITE
233  ),
234  mesh_,
236  ),
237  qr_
238  (
239  IOobject
240  (
241  "qr",
242  mesh_.time().timeName(),
243  mesh_,
244  IOobject::READ_IF_PRESENT,
245  IOobject::AUTO_WRITE
246  ),
247  mesh_,
249  ),
250  qem_
251  (
252  IOobject
253  (
254  "qem",
255  mesh_.time().timeName(),
256  mesh_,
257  IOobject::READ_IF_PRESENT,
258  IOobject::AUTO_WRITE
259  ),
260  mesh_,
262  ),
263  qin_
264  (
265  IOobject
266  (
267  "qin",
268  mesh_.time().timeName(),
269  mesh_,
270  IOobject::READ_IF_PRESENT,
271  IOobject::AUTO_WRITE
272  ),
273  mesh_,
275  ),
276  a_
277  (
278  IOobject
279  (
280  "a",
281  mesh_.time().timeName(),
282  mesh_,
283  IOobject::NO_READ,
284  IOobject::AUTO_WRITE
285  ),
286  mesh_,
288  ),
289  nTheta_(coeffs_.lookup<label>("nTheta")),
290  nPhi_(coeffs_.lookup<label>("nPhi")),
291  nRay_(0),
292  nLambda_(absorptionEmission_->nBands()),
293  aLambda_(nLambda_),
294  blackBody_(nLambda_, T),
295  IRay_(0),
296  tolerance_
297  (
298  coeffs_.lookupOrDefaultBackwardsCompatible<scalar>
299  (
300  {"tolerance", "convergence"},
301  0
302  )
303  ),
304  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
305  omegaMax_(0)
306 {
307  initialise();
308 }
309 
310 
312 (
313  const dictionary& dict,
314  const volScalarField& T
315 )
316 :
317  radiationModel(typeName, dict, T),
318  G_
319  (
320  IOobject
321  (
322  "G",
323  mesh_.time().timeName(),
324  mesh_,
327  ),
328  mesh_,
330  ),
331  qr_
332  (
333  IOobject
334  (
335  "qr",
336  mesh_.time().timeName(),
337  mesh_,
340  ),
341  mesh_,
343  ),
344  qem_
345  (
346  IOobject
347  (
348  "qem",
349  mesh_.time().timeName(),
350  mesh_,
353  ),
354  mesh_,
356  ),
357  qin_
358  (
359  IOobject
360  (
361  "qin",
362  mesh_.time().timeName(),
363  mesh_,
366  ),
367  mesh_,
369  ),
370  a_
371  (
372  IOobject
373  (
374  "a",
375  mesh_.time().timeName(),
376  mesh_,
379  ),
380  mesh_,
382  ),
383  nTheta_(coeffs_.lookup<label>("nTheta")),
384  nPhi_(coeffs_.lookup<label>("nPhi")),
385  nRay_(0),
386  nLambda_(absorptionEmission_->nBands()),
387  aLambda_(nLambda_),
388  blackBody_(nLambda_, T),
389  IRay_(0),
390  tolerance_
391  (
393  (
394  {"tolerance", "convergence"},
395  0
396  )
397  ),
398  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
399  omegaMax_(0)
400 {
401  initialise();
402 }
403 
404 
405 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
406 
408 {}
409 
410 
411 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
412 
414 {
415  if (radiationModel::read())
416  {
417  // Only reading solution parameters - not changing ray geometry
418 
419  // For backward-compatibility
420  coeffs_.readIfPresent("convergence", tolerance_);
421  coeffs_.readIfPresent("tolerance", tolerance_);
422  coeffs_.readIfPresent("maxIter", maxIter_);
423 
424  return true;
425  }
426  else
427  {
428  return false;
429  }
430 }
431 
432 
434 {
435  absorptionEmission_->correct(a_, aLambda_);
436 
437  updateBlackBodyEmission();
438 
439  // Set rays converged false
440  List<bool> rayIdConv(nRay_, false);
441 
442  scalar maxResidual = 0;
443  label radIter = 0;
444  do
445  {
446  Info<< "Radiation solver iter: " << radIter << endl;
447 
448  radIter++;
449  maxResidual = 0;
450  forAll(IRay_, rayI)
451  {
452  if (!rayIdConv[rayI])
453  {
454  scalar maxBandResidual = IRay_[rayI].correct();
455  maxResidual = max(maxBandResidual, maxResidual);
456 
457  if (maxBandResidual < tolerance_)
458  {
459  rayIdConv[rayI] = true;
460  }
461  }
462  }
463 
464  } while (maxResidual > tolerance_ && radIter < maxIter_);
465 
466  updateG();
467 }
468 
469 
471 {
472  // Construct using contribution from first frequency band
474  (
476  (
477  "Rp",
478  4
480  *(aLambda_[0] - absorptionEmission_->aDisp(0)())
481  *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(0))
482  )
483  );
484 
485  volScalarField& Rp=tRp.ref();
486 
487  // Add contributions over remaining frequency bands
488  for (label j=1; j < nLambda_; j++)
489  {
490  Rp +=
491  (
492  4
494  *(aLambda_[j] - absorptionEmission_->aDisp(j)())
495  *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(j))
496  );
497  }
498 
499  return tRp;
500 }
501 
502 
505 {
507  (
509  (
510  IOobject
511  (
512  "Ru",
513  mesh_.time().timeName(),
514  mesh_,
517  false
518  ),
519  mesh_,
520  dimensionedScalar(dimensionSet(1, -1, -3, 0, 0), 0)
521  )
522  );
523 
525 
526  // Sum contributions over all frequency bands
527  for (label j=0; j < nLambda_; j++)
528  {
529  // Compute total incident radiation within frequency band
531  (
532  IRay_[0].ILambda(j)()*IRay_[0].omega()
533  );
534 
535  for (label rayI=1; rayI < nRay_; rayI++)
536  {
537  Gj.ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
538  }
539 
540  Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
541  - absorptionEmission_->ECont(j)()();
542  }
543 
544  return tRu;
545 }
546 
547 
548 void Foam::radiationModels::fvDOM::updateBlackBodyEmission()
549 {
550  for (label j=0; j < nLambda_; j++)
551  {
552  blackBody_.correct(j, absorptionEmission_->bands(j));
553  }
554 }
555 
556 
558 {
559  G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0);
560  qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0);
563 
564  forAll(IRay_, rayI)
565  {
566  IRay_[rayI].addIntensity();
567  G_ += IRay_[rayI].I()*IRay_[rayI].omega();
568  qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
569  qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
570  qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
571  }
572 }
573 
574 
576 (
577  const word& name,
578  label& rayId,
579  label& lambdaId
580 ) const
581 {
582  // Assume name is in the form: <name>_<rayId>_<lambdaId>
583  const size_type i1 = name.find_first_of("_");
584  const size_type i2 = name.find_last_of("_");
585 
586  rayId = readLabel(IStringStream(name.substr(i1 + 1, i2 - 1))());
587  lambdaId = readLabel(IStringStream(name.substr(i2 + 1, name.size() - 1))());
588 }
589 
590 
591 // ************************************************************************* //
Collection of constants.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void calculate()
Solve radiation equation(s)
Definition: fvDOM.C:433
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
tmp< Foam::volScalarField > deltaLambdaT(const volScalarField &T, const Vector2D< scalar > &band) const
Proportion of total energy at T from lambda1 to lambda2.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual bool read()=0
Read radiationProperties dictionary.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:636
const dimensionSet dimless
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
Macros for easy insertion into run-time selection tables.
const dimensionSet dimLength
virtual tmp< volScalarField::Internal > Ru() const
Source term component (constant)
Definition: fvDOM.C:504
const volScalarField & T_
Reference to the temperature field.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
fvDOM(const volScalarField &T)
Construct from components.
Definition: fvDOM.C:221
const dimensionSet dimTime
void correct(const label lambdaI, const Vector2D< scalar > &band)
T lookupOrDefaultBackwardsCompatible(const wordList &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T, trying a list of keywords in sequence.
Dimension set for the base types.
Definition: dimensionSet.H:120
stressControl lookup("compactNormalStress") >> compactNormalStress
mathematical constants.
A class for handling words, derived from string.
Definition: word.H:59
Top level model for radiation modelling.
virtual const fileName & name() const
Return the name of the stream.
Definition: IOstream.H:294
bool read()
Read radiation properties dictionary.
Definition: fvDOM.C:413
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
word timeName
Definition: getTimeIndex.H:3
label readLabel(Istream &is)
Definition: label.H:64
void setRayIdLambdaId(const word &name, label &rayId, label &lambdaId) const
Set the rayId and lambdaId from by decomposing an intensity.
Definition: fvDOM.C:576
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:235
const dimensionSet dimMass
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
dimensionedScalar pow3(const dimensionedScalar &ds)
label patchi
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label size_type
The type that can represent the size of a DLList.
Definition: UILList.H:177
Input from memory buffer stream.
Definition: IStringStream.H:49
radiationModel(const volScalarField &T)
Null constructor.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
messageStream Info
const fvMesh & mesh_
Reference to the mesh database.
label n
virtual ~fvDOM()
Destructor.
Definition: fvDOM.C:407
A class for managing temporary objects.
Definition: PtrList.H:53
const scalar piByTwo(0.5 *pi)
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:228
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
#define addToRadiationRunTimeSelectionTables(model)
autoPtr< radiationModels::absorptionEmissionModel > absorptionEmission_
Absorption/emission model.
virtual tmp< volScalarField > Rp() const
Source term component (for power of T^4)
Definition: fvDOM.C:470
dictionary coeffs_
Radiation model dictionary.
void updateG()
Update G and calculate total heat flux on boundary.
Definition: fvDOM.C:557
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844