All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2019 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"
32 
33 using namespace Foam::constant;
34 using namespace Foam::constant::mathematical;
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 namespace radiationModels
41 {
42  defineTypeNameAndDebug(fvDOM, 0);
44 }
45 }
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 void Foam::radiationModels::fvDOM::initialise()
51 {
52  // 3D
53  if (mesh_.nSolutionD() == 3)
54  {
55  nRay_ = 4*nPhi_*nTheta_;
56  IRay_.setSize(nRay_);
57  const scalar deltaPhi = pi/(2.0*nPhi_);
58  const scalar deltaTheta = pi/nTheta_;
59  label i = 0;
60  for (label n = 1; n <= nTheta_; n++)
61  {
62  for (label m = 1; m <= 4*nPhi_; m++)
63  {
64  const scalar thetai = (2*n - 1)*deltaTheta/2.0;
65  const scalar phii = (2*m - 1)*deltaPhi/2.0;
66  IRay_.set
67  (
68  i,
69  new radiativeIntensityRay
70  (
71  *this,
72  mesh_,
73  phii,
74  thetai,
75  deltaPhi,
76  deltaTheta,
77  nLambda_,
78  absorptionEmission_,
79  blackBody_,
80  i
81  )
82  );
83  i++;
84  }
85  }
86  }
87  // 2D
88  else if (mesh_.nSolutionD() == 2)
89  {
90  const scalar thetai = piByTwo;
91  const scalar deltaTheta = pi;
92  nRay_ = 4*nPhi_;
93  IRay_.setSize(nRay_);
94  const scalar deltaPhi = pi/(2.0*nPhi_);
95  label i = 0;
96  for (label m = 1; m <= 4*nPhi_; m++)
97  {
98  const scalar phii = (2*m - 1)*deltaPhi/2.0;
99  IRay_.set
100  (
101  i,
102  new radiativeIntensityRay
103  (
104  *this,
105  mesh_,
106  phii,
107  thetai,
108  deltaPhi,
109  deltaTheta,
110  nLambda_,
111  absorptionEmission_,
112  blackBody_,
113  i
114  )
115  );
116  i++;
117  }
118  }
119  // 1D
120  else
121  {
122  const scalar thetai = piByTwo;
123  const scalar deltaTheta = pi;
124  nRay_ = 2;
125  IRay_.setSize(nRay_);
126  const scalar deltaPhi = pi;
127  label i = 0;
128  for (label m = 1; m <= 2; m++)
129  {
130  const scalar phii = (2*m - 1)*deltaPhi/2.0;
131  IRay_.set
132  (
133  i,
134  new radiativeIntensityRay
135  (
136  *this,
137  mesh_,
138  phii,
139  thetai,
140  deltaPhi,
141  deltaTheta,
142  nLambda_,
143  absorptionEmission_,
144  blackBody_,
145  i
146  )
147  );
148  i++;
149  }
150  }
151 
152 
153  // Construct absorption field for each wavelength
154  forAll(aLambda_, lambdaI)
155  {
156  aLambda_.set
157  (
158  lambdaI,
159  new volScalarField
160  (
161  IOobject
162  (
163  "aLambda_" + Foam::name(lambdaI) ,
164  mesh_.time().timeName(),
165  mesh_,
168  ),
169  a_
170  )
171  );
172  }
173 
174 
175  // Calculate the maximum solid angle
176  forAll(IRay_, rayId)
177  {
178  if (omegaMax_ < IRay_[rayId].omega())
179  {
180  omegaMax_ = IRay_[rayId].omega();
181  }
182  }
183 
184 
185  Info<< typeName << ": Created " << IRay_.size() << " rays with average "
186  << "directions (dAve) and solid angles (omega)" << endl;
187  Info<< incrIndent;
188  forAll(IRay_, rayId)
189  {
190  Info<< indent
191  << "Ray " << IRay_[rayId].I().name() << ": "
192  << "dAve = " << IRay_[rayId].dAve() << ", "
193  << "omega = " << IRay_[rayId].omega() << endl;
194  }
195  Info<< decrIndent << endl;
196 }
197 
198 
199 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
200 
202 :
203  radiationModel(typeName, T),
204  G_
205  (
206  IOobject
207  (
208  "G",
209  mesh_.time().timeName(),
210  mesh_,
211  IOobject::NO_READ,
212  IOobject::AUTO_WRITE
213  ),
214  mesh_,
216  ),
217  qr_
218  (
219  IOobject
220  (
221  "qr",
222  mesh_.time().timeName(),
223  mesh_,
224  IOobject::READ_IF_PRESENT,
225  IOobject::AUTO_WRITE
226  ),
227  mesh_,
229  ),
230  qem_
231  (
232  IOobject
233  (
234  "qem",
235  mesh_.time().timeName(),
236  mesh_,
237  IOobject::READ_IF_PRESENT,
238  IOobject::AUTO_WRITE
239  ),
240  mesh_,
242  ),
243  qin_
244  (
245  IOobject
246  (
247  "qin",
248  mesh_.time().timeName(),
249  mesh_,
250  IOobject::READ_IF_PRESENT,
251  IOobject::AUTO_WRITE
252  ),
253  mesh_,
255  ),
256  a_
257  (
258  IOobject
259  (
260  "a",
261  mesh_.time().timeName(),
262  mesh_,
263  IOobject::NO_READ,
264  IOobject::AUTO_WRITE
265  ),
266  mesh_,
268  ),
269  nTheta_(coeffs_.lookup<label>("nTheta")),
270  nPhi_(coeffs_.lookup<label>("nPhi")),
271  nRay_(0),
272  nLambda_(absorptionEmission_->nBands()),
273  aLambda_(nLambda_),
274  blackBody_(nLambda_, T),
275  IRay_(0),
276  tolerance_
277  (
278  coeffs_.found("convergence")
279  ? coeffs_.lookup<scalar>("convergence")
280  : coeffs_.lookupOrDefault<scalar>("tolerance", 0)
281  ),
282  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
283  omegaMax_(0)
284 {
285  initialise();
286 }
287 
288 
290 (
291  const dictionary& dict,
292  const volScalarField& T
293 )
294 :
295  radiationModel(typeName, dict, T),
296  G_
297  (
298  IOobject
299  (
300  "G",
301  mesh_.time().timeName(),
302  mesh_,
305  ),
306  mesh_,
308  ),
309  qr_
310  (
311  IOobject
312  (
313  "qr",
314  mesh_.time().timeName(),
315  mesh_,
318  ),
319  mesh_,
321  ),
322  qem_
323  (
324  IOobject
325  (
326  "qem",
327  mesh_.time().timeName(),
328  mesh_,
331  ),
332  mesh_,
334  ),
335  qin_
336  (
337  IOobject
338  (
339  "qin",
340  mesh_.time().timeName(),
341  mesh_,
344  ),
345  mesh_,
347  ),
348  a_
349  (
350  IOobject
351  (
352  "a",
353  mesh_.time().timeName(),
354  mesh_,
357  ),
358  mesh_,
360  ),
361  nTheta_(coeffs_.lookup<label>("nTheta")),
362  nPhi_(coeffs_.lookup<label>("nPhi")),
363  nRay_(0),
364  nLambda_(absorptionEmission_->nBands()),
365  aLambda_(nLambda_),
366  blackBody_(nLambda_, T),
367  IRay_(0),
368  tolerance_
369  (
370  coeffs_.found("convergence")
371  ? coeffs_.lookup<scalar>("convergence")
372  : coeffs_.lookupOrDefault<scalar>("tolerance", 0)
373  ),
374  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
375  omegaMax_(0)
376 {
377  initialise();
378 }
379 
380 
381 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
382 
384 {}
385 
386 
387 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
388 
390 {
391  if (radiationModel::read())
392  {
393  // Only reading solution parameters - not changing ray geometry
394 
395  // For backward-compatibility
396  coeffs_.readIfPresent("convergence", tolerance_);
397  coeffs_.readIfPresent("tolerance", tolerance_);
398  coeffs_.readIfPresent("maxIter", maxIter_);
399 
400  return true;
401  }
402  else
403  {
404  return false;
405  }
406 }
407 
408 
410 {
411  absorptionEmission_->correct(a_, aLambda_);
412 
413  updateBlackBodyEmission();
414 
415  // Set rays converged false
416  List<bool> rayIdConv(nRay_, false);
417 
418  scalar maxResidual = 0;
419  label radIter = 0;
420  do
421  {
422  Info<< "Radiation solver iter: " << radIter << endl;
423 
424  radIter++;
425  maxResidual = 0;
426  forAll(IRay_, rayI)
427  {
428  if (!rayIdConv[rayI])
429  {
430  scalar maxBandResidual = IRay_[rayI].correct();
431  maxResidual = max(maxBandResidual, maxResidual);
432 
433  if (maxBandResidual < tolerance_)
434  {
435  rayIdConv[rayI] = true;
436  }
437  }
438  }
439 
440  } while (maxResidual > tolerance_ && radIter < maxIter_);
441 
442  updateG();
443 }
444 
445 
447 {
448  // Construct using contribution from first frequency band
450  (
452  (
453  "Rp",
454  4
456  *(aLambda_[0] - absorptionEmission_->aDisp(0)())
457  *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(0))
458  )
459  );
460 
461  volScalarField& Rp=tRp.ref();
462 
463  // Add contributions over remaining frequency bands
464  for (label j=1; j < nLambda_; j++)
465  {
466  Rp +=
467  (
468  4
470  *(aLambda_[j] - absorptionEmission_->aDisp(j)())
471  *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(j))
472  );
473  }
474 
475  return tRp;
476 }
477 
478 
481 {
483  (
485  (
486  IOobject
487  (
488  "Ru",
489  mesh_.time().timeName(),
490  mesh_,
493  false
494  ),
495  mesh_,
496  dimensionedScalar(dimensionSet(1, -1, -3, 0, 0), 0)
497  )
498  );
499 
501 
502  // Sum contributions over all frequency bands
503  for (label j=0; j < nLambda_; j++)
504  {
505  // Compute total incident radiation within frequency band
507  (
508  IRay_[0].ILambda(j)()*IRay_[0].omega()
509  );
510 
511  for (label rayI=1; rayI < nRay_; rayI++)
512  {
513  Gj.ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
514  }
515 
516  Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
517  - absorptionEmission_->ECont(j)()();
518  }
519 
520  return tRu;
521 }
522 
523 
524 void Foam::radiationModels::fvDOM::updateBlackBodyEmission()
525 {
526  for (label j=0; j < nLambda_; j++)
527  {
528  blackBody_.correct(j, absorptionEmission_->bands(j));
529  }
530 }
531 
532 
534 {
535  G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0);
536  qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0);
539 
540  forAll(IRay_, rayI)
541  {
542  IRay_[rayI].addIntensity();
543  G_ += IRay_[rayI].I()*IRay_[rayI].omega();
544  qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
545  qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
546  qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
547  }
548 }
549 
550 
552 (
553  const word& name,
554  label& rayId,
555  label& lambdaId
556 ) const
557 {
558  // Assume name is in the form: <name>_<rayId>_<lambdaId>
559  const size_type i1 = name.find_first_of("_");
560  const size_type i2 = name.find_last_of("_");
561 
562  rayId = readLabel(IStringStream(name.substr(i1 + 1, i2 - 1))());
563  lambdaId = readLabel(IStringStream(name.substr(i2 + 1, name.size() - 1))());
564 }
565 
566 
567 // ************************************************************************* //
Collection of constants.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:667
#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:409
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
tmp< Foam::volScalarField > deltaLambdaT(const volScalarField &T, const Vector2D< scalar > &band) const
Proportion of total energy at T from lambda1 to lambda2.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
const dimensionedScalar & sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
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:622
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
Macros for easy insertion into run-time selection tables.
virtual tmp< volScalarField::Internal > Ru() const
Source term component (constant)
Definition: fvDOM.C:480
const volScalarField & T_
Reference to the temperature field.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
fvDOM(const volScalarField &T)
Construct from components.
Definition: fvDOM.C:201
void correct(const label lambdaI, const Vector2D< scalar > &band)
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:297
bool read()
Read radiation properties dictionary.
Definition: fvDOM.C:389
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:552
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:235
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
dimensionedScalar pow3(const dimensionedScalar &ds)
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
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.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
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.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
messageStream Info
const fvMesh & mesh_
Reference to the mesh database.
label n
virtual ~fvDOM()
Destructor.
Definition: fvDOM.C:383
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
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.
bool found
virtual tmp< volScalarField > Rp() const
Source term component (for power of T^4)
Definition: fvDOM.C:446
dictionary coeffs_
Radiation model dictionary.
void updateG()
Update G and calculate total heat flux on boundary.
Definition: fvDOM.C:533
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:812