fvDOM.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 radiation
41  {
42  defineTypeNameAndDebug(fvDOM, 0);
44  }
45 }
46 
47 
48 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
49 
50 void Foam::radiation::fvDOM::initialise()
51 {
52  // 3D
53  if (mesh_.nSolutionD() == 3)
54  {
55  nRay_ = 4*nPhi_*nTheta_;
56  IRay_.setSize(nRay_);
57  scalar deltaPhi = pi/(2.0*nPhi_);
58  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  scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
65  scalar phii = (2.0*m - 1.0)*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  // Currently 2D solution is limited to the x-y plane
91  if (mesh_.solutionD()[vector::Z] != -1)
92  {
94  << "Currently 2D solution is limited to the x-y plane"
95  << exit(FatalError);
96  }
97 
98  scalar thetai = piByTwo;
99  scalar deltaTheta = pi;
100  nRay_ = 4*nPhi_;
101  IRay_.setSize(nRay_);
102  scalar deltaPhi = pi/(2.0*nPhi_);
103  label i = 0;
104  for (label m = 1; m <= 4*nPhi_; m++)
105  {
106  scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
107  IRay_.set
108  (
109  i,
110  new radiativeIntensityRay
111  (
112  *this,
113  mesh_,
114  phii,
115  thetai,
116  deltaPhi,
117  deltaTheta,
118  nLambda_,
119  absorptionEmission_,
120  blackBody_,
121  i
122  )
123  );
124  i++;
125  }
126  }
127  // 1D
128  else
129  {
130  // Currently 1D solution is limited to the x-direction
131  if (mesh_.solutionD()[vector::X] != 1)
132  {
134  << "Currently 1D solution is limited to the x-direction"
135  << exit(FatalError);
136  }
137 
138  scalar thetai = piByTwo;
139  scalar deltaTheta = pi;
140  nRay_ = 2;
141  IRay_.setSize(nRay_);
142  scalar deltaPhi = pi;
143  label i = 0;
144  for (label m = 1; m <= 2; m++)
145  {
146  scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
147  IRay_.set
148  (
149  i,
150  new radiativeIntensityRay
151  (
152  *this,
153  mesh_,
154  phii,
155  thetai,
156  deltaPhi,
157  deltaTheta,
158  nLambda_,
159  absorptionEmission_,
160  blackBody_,
161  i
162  )
163  );
164  i++;
165  }
166  }
167 
168 
169  // Construct absorption field for each wavelength
170  forAll(aLambda_, lambdaI)
171  {
172  aLambda_.set
173  (
174  lambdaI,
175  new volScalarField
176  (
177  IOobject
178  (
179  "aLambda_" + Foam::name(lambdaI) ,
180  mesh_.time().timeName(),
181  mesh_,
184  ),
185  a_
186  )
187  );
188  }
189 
190  Info<< "fvDOM : Allocated " << IRay_.size()
191  << " rays with average orientation:" << nl;
192 
193  if (cacheDiv_)
194  {
195  Info<< "Caching div fvMatrix..."<< endl;
196  for (label lambdaI = 0; lambdaI < nLambda_; lambdaI++)
197  {
198  fvRayDiv_[lambdaI].setSize(nRay_);
199 
200  forAll(IRay_, rayId)
201  {
202  const surfaceScalarField Ji(IRay_[rayId].dAve() & mesh_.Sf());
203  const volScalarField& iRayLambdaI =
204  IRay_[rayId].ILambda(lambdaI);
205 
206  fvRayDiv_[lambdaI].set
207  (
208  rayId,
209  new fvScalarMatrix
210  (
211  fvm::div(Ji, iRayLambdaI, "div(Ji,Ii_h)")
212  )
213  );
214  }
215  }
216  }
217 
218  forAll(IRay_, rayId)
219  {
220  if (omegaMax_ < IRay_[rayId].omega())
221  {
222  omegaMax_ = IRay_[rayId].omega();
223  }
224  Info<< '\t' << IRay_[rayId].I().name() << " : " << "omega : "
225  << '\t' << IRay_[rayId].omega() << nl;
226  }
227 
228  Info<< endl;
229 }
230 
231 
232 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
233 
234 Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
235 :
236  radiationModel(typeName, T),
237  G_
238  (
239  IOobject
240  (
241  "G",
242  mesh_.time().timeName(),
243  mesh_,
244  IOobject::NO_READ,
245  IOobject::AUTO_WRITE
246  ),
247  mesh_,
249  ),
250  Qr_
251  (
252  IOobject
253  (
254  "Qr",
255  mesh_.time().timeName(),
256  mesh_,
257  IOobject::READ_IF_PRESENT,
258  IOobject::AUTO_WRITE
259  ),
260  mesh_,
262  ),
263  Qem_
264  (
265  IOobject
266  (
267  "Qem",
268  mesh_.time().timeName(),
269  mesh_,
270  IOobject::NO_READ,
271  IOobject::NO_WRITE
272  ),
273  mesh_,
274  dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
275  ),
276  Qin_
277  (
278  IOobject
279  (
280  "Qin",
281  mesh_.time().timeName(),
282  mesh_,
283  IOobject::READ_IF_PRESENT,
284  IOobject::AUTO_WRITE
285  ),
286  mesh_,
287  dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
288  ),
289  a_
290  (
291  IOobject
292  (
293  "a",
294  mesh_.time().timeName(),
295  mesh_,
296  IOobject::NO_READ,
297  IOobject::AUTO_WRITE
298  ),
299  mesh_,
301  ),
302  nTheta_(readLabel(coeffs_.lookup("nTheta"))),
303  nPhi_(readLabel(coeffs_.lookup("nPhi"))),
304  nRay_(0),
305  nLambda_(absorptionEmission_->nBands()),
306  aLambda_(nLambda_),
307  blackBody_(nLambda_, T),
308  IRay_(0),
309  convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
310  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
311  fvRayDiv_(nLambda_),
312  cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
313  omegaMax_(0)
314 {
315  initialise();
316 }
317 
318 
319 Foam::radiation::fvDOM::fvDOM
320 (
321  const dictionary& dict,
322  const volScalarField& T
323 )
324 :
325  radiationModel(typeName, dict, T),
326  G_
327  (
328  IOobject
329  (
330  "G",
331  mesh_.time().timeName(),
332  mesh_,
335  ),
336  mesh_,
338  ),
339  Qr_
340  (
341  IOobject
342  (
343  "Qr",
344  mesh_.time().timeName(),
345  mesh_,
348  ),
349  mesh_,
351  ),
352  Qem_
353  (
354  IOobject
355  (
356  "Qem",
357  mesh_.time().timeName(),
358  mesh_,
361  ),
362  mesh_,
363  dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
364  ),
365  Qin_
366  (
367  IOobject
368  (
369  "Qin",
370  mesh_.time().timeName(),
371  mesh_,
374  ),
375  mesh_,
376  dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
377  ),
378  a_
379  (
380  IOobject
381  (
382  "a",
383  mesh_.time().timeName(),
384  mesh_,
387  ),
388  mesh_,
390  ),
391  nTheta_(readLabel(coeffs_.lookup("nTheta"))),
392  nPhi_(readLabel(coeffs_.lookup("nPhi"))),
393  nRay_(0),
394  nLambda_(absorptionEmission_->nBands()),
395  aLambda_(nLambda_),
396  blackBody_(nLambda_, T),
397  IRay_(0),
398  convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
399  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
400  fvRayDiv_(nLambda_),
401  cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
402  omegaMax_(0)
403 {
404  initialise();
405 }
406 
407 
408 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
409 
411 {}
412 
413 
414 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
415 
417 {
418  if (radiationModel::read())
419  {
420  // Only reading solution parameters - not changing ray geometry
421  coeffs_.readIfPresent("convergence", convergence_);
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 convergence false
440  List<bool> rayIdConv(nRay_, false);
441 
442  scalar maxResidual = 0.0;
443  label radIter = 0;
444  do
445  {
446  Info<< "Radiation solver iter: " << radIter << endl;
447 
448  radIter++;
449  maxResidual = 0.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 < convergence_)
458  {
459  rayIdConv[rayI] = true;
460  }
461  }
462  }
463 
464  } while (maxResidual > convergence_ && radIter < maxIter_);
465 
466  updateG();
467 }
468 
469 
471 {
472  return tmp<volScalarField>
473  (
474  new volScalarField
475  (
476  IOobject
477  (
478  "Rp",
479  mesh_.time().timeName(),
480  mesh_,
483  false
484  ),
485  // Only include continuous phase emission
487  )
488  );
489 }
490 
491 
494 {
495 
497  G_();
498 
500  absorptionEmission_->ECont()()();
501 
502  // Only include continuous phase absorption
504  absorptionEmission_->aCont()()();
505 
506  return a*G - E;
507 }
508 
509 
510 void Foam::radiation::fvDOM::updateBlackBodyEmission()
511 {
512  for (label j=0; j < nLambda_; j++)
513  {
514  blackBody_.correct(j, absorptionEmission_->bands(j));
515  }
516 }
517 
518 
520 {
521  G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
522  Qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
523  Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
524  Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
525 
526  forAll(IRay_, rayI)
527  {
528  IRay_[rayI].addIntensity();
529  G_ += IRay_[rayI].I()*IRay_[rayI].omega();
530  Qr_.boundaryFieldRef() += IRay_[rayI].Qr().boundaryField();
531  Qem_.boundaryFieldRef() += IRay_[rayI].Qem().boundaryField();
532  Qin_.boundaryFieldRef() += IRay_[rayI].Qin().boundaryField();
533  }
534 }
535 
536 
538 (
539  const word& name,
540  label& rayId,
541  label& lambdaId
542 ) const
543 {
544  // assuming name is in the form: CHARS_rayId_lambdaId
545  size_type i1 = name.find_first_of("_");
546  size_type i2 = name.find_last_of("_");
547 
548  rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
549  lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
550 }
551 
552 
553 // ************************************************************************* //
Collection of constants.
dictionary coeffs_
Radiation model dictionary.
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
void correct(const label lambdaI, const Vector2D< scalar > &band)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
const volScalarField & G() const
Const access to incident radiation field.
Definition: fvDOM.H:84
void setRayIdLambdaId(const word &name, label &rayId, label &lambdaId) const
Set the rayId and lambdaId from by decomposing an intensity.
Definition: fvDOM.C:538
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:715
Macros for easy insertion into run-time selection tables.
virtual tmp< DimensionedField< scalar, volMesh > > Ru() const
Source term component (constant)
Definition: fvDOM.C:493
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
stressControl lookup("compactNormalStress") >> compactNormalStress
virtual bool read()=0
Read radiationProperties dictionary.
mathematical constants.
A class for handling words, derived from string.
Definition: word.H:59
void updateG()
Update G and calculate total heat flux on boundary.
Definition: fvDOM.C:519
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
const fvMesh & mesh_
Reference to the mesh database.
Top level model for radiation modelling.
word timeName
Definition: getTimeIndex.H:3
label readLabel(Istream &is)
Definition: label.H:64
bool read()
Read radiation properties dictionary.
Definition: fvDOM.C:416
static const char nl
Definition: Ostream.H:262
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
const volScalarField & a() const
Const access to total absorption coefficient.
Definition: fvDOM.H:69
defineTypeNameAndDebug(combustionModel, 0)
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
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:170
virtual ~fvDOM()
Destructor.
Definition: fvDOM.C:410
Input from memory buffer stream.
Definition: IStringStream.H:49
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
void calculate()
Solve radiation equation(s)
Definition: fvDOM.C:433
messageStream Info
label n
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
A class for managing temporary objects.
Definition: PtrList.H:54
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
const scalar piByTwo(0.5 *pi)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define addToRadiationRunTimeSelectionTables(model)
autoPtr< absorptionEmissionModel > absorptionEmission_
Absorption/emission model.
autoPtr< radiation::radiationModel > radiation(radiation::radiationModel::New(T))
Namespace for OpenFOAM.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
virtual tmp< volScalarField > Rp() const
Source term component (for power of T^4)
Definition: fvDOM.C:470
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451