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-2017 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  forAll(IRay_, rayId)
194  {
195  if (omegaMax_ < IRay_[rayId].omega())
196  {
197  omegaMax_ = IRay_[rayId].omega();
198  }
199  Info<< '\t' << IRay_[rayId].I().name() << " : " << "omega : "
200  << '\t' << IRay_[rayId].omega() << nl;
201  }
202 
203  Info<< endl;
204 }
205 
206 
207 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
208 
209 Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
210 :
211  radiationModel(typeName, T),
212  G_
213  (
214  IOobject
215  (
216  "G",
217  mesh_.time().timeName(),
218  mesh_,
219  IOobject::NO_READ,
220  IOobject::AUTO_WRITE
221  ),
222  mesh_,
224  ),
225  qr_
226  (
227  IOobject
228  (
229  "qr",
230  mesh_.time().timeName(),
231  mesh_,
232  IOobject::READ_IF_PRESENT,
233  IOobject::AUTO_WRITE
234  ),
235  mesh_,
237  ),
238  qem_
239  (
240  IOobject
241  (
242  "qem",
243  mesh_.time().timeName(),
244  mesh_,
245  IOobject::NO_READ,
246  IOobject::NO_WRITE
247  ),
248  mesh_,
249  dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
250  ),
251  qin_
252  (
253  IOobject
254  (
255  "qin",
256  mesh_.time().timeName(),
257  mesh_,
258  IOobject::READ_IF_PRESENT,
259  IOobject::AUTO_WRITE
260  ),
261  mesh_,
262  dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
263  ),
264  a_
265  (
266  IOobject
267  (
268  "a",
269  mesh_.time().timeName(),
270  mesh_,
271  IOobject::NO_READ,
272  IOobject::AUTO_WRITE
273  ),
274  mesh_,
276  ),
277  nTheta_(readLabel(coeffs_.lookup("nTheta"))),
278  nPhi_(readLabel(coeffs_.lookup("nPhi"))),
279  nRay_(0),
280  nLambda_(absorptionEmission_->nBands()),
281  aLambda_(nLambda_),
282  blackBody_(nLambda_, T),
283  IRay_(0),
284  tolerance_
285  (
286  coeffs_.found("convergence")
287  ? readScalar(coeffs_.lookup("convergence"))
288  : coeffs_.lookupOrDefault<scalar>("tolerance", 0.0)
289  ),
290  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
291  omegaMax_(0)
292 {
293  initialise();
294 }
295 
296 
297 Foam::radiation::fvDOM::fvDOM
298 (
299  const dictionary& dict,
300  const volScalarField& T
301 )
302 :
303  radiationModel(typeName, dict, T),
304  G_
305  (
306  IOobject
307  (
308  "G",
309  mesh_.time().timeName(),
310  mesh_,
313  ),
314  mesh_,
316  ),
317  qr_
318  (
319  IOobject
320  (
321  "qr",
322  mesh_.time().timeName(),
323  mesh_,
326  ),
327  mesh_,
329  ),
330  qem_
331  (
332  IOobject
333  (
334  "qem",
335  mesh_.time().timeName(),
336  mesh_,
339  ),
340  mesh_,
341  dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
342  ),
343  qin_
344  (
345  IOobject
346  (
347  "qin",
348  mesh_.time().timeName(),
349  mesh_,
352  ),
353  mesh_,
354  dimensionedScalar("qin", dimMass/pow3(dimTime), 0.0)
355  ),
356  a_
357  (
358  IOobject
359  (
360  "a",
361  mesh_.time().timeName(),
362  mesh_,
365  ),
366  mesh_,
368  ),
369  nTheta_(readLabel(coeffs_.lookup("nTheta"))),
370  nPhi_(readLabel(coeffs_.lookup("nPhi"))),
371  nRay_(0),
372  nLambda_(absorptionEmission_->nBands()),
373  aLambda_(nLambda_),
374  blackBody_(nLambda_, T),
375  IRay_(0),
376  tolerance_
377  (
378  coeffs_.found("convergence")
379  ? readScalar(coeffs_.lookup("convergence"))
380  : coeffs_.lookupOrDefault<scalar>("tolerance", 0.0)
381  ),
382  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
383  omegaMax_(0)
384 {
385  initialise();
386 }
387 
388 
389 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
390 
392 {}
393 
394 
395 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
396 
398 {
399  if (radiationModel::read())
400  {
401  // Only reading solution parameters - not changing ray geometry
402 
403  // For backward-compatibility
404  coeffs_.readIfPresent("convergence", tolerance_);
405 
406  coeffs_.readIfPresent("tolerance", tolerance_);
407 
408  coeffs_.readIfPresent("maxIter", maxIter_);
409 
410  return true;
411  }
412  else
413  {
414  return false;
415  }
416 }
417 
418 
420 {
421  absorptionEmission_->correct(a_, aLambda_);
422 
423  updateBlackBodyEmission();
424 
425  // Set rays converged false
426  List<bool> rayIdConv(nRay_, false);
427 
428  scalar maxResidual = 0.0;
429  label radIter = 0;
430  do
431  {
432  Info<< "Radiation solver iter: " << radIter << endl;
433 
434  radIter++;
435  maxResidual = 0.0;
436  forAll(IRay_, rayI)
437  {
438  if (!rayIdConv[rayI])
439  {
440  scalar maxBandResidual = IRay_[rayI].correct();
441  maxResidual = max(maxBandResidual, maxResidual);
442 
443  if (maxBandResidual < tolerance_)
444  {
445  rayIdConv[rayI] = true;
446  }
447  }
448  }
449 
450  } while (maxResidual > tolerance_ && radIter < maxIter_);
451 
452  updateG();
453 }
454 
455 
457 {
458  return tmp<volScalarField>
459  (
460  new volScalarField
461  (
462  IOobject
463  (
464  "Rp",
465  mesh_.time().timeName(),
466  mesh_,
469  false
470  ),
471  // Only include continuous phase emission
473  )
474  );
475 }
476 
477 
480 {
481 
482  const volScalarField::Internal& G =
483  G_();
484 
485  const volScalarField::Internal E =
486  absorptionEmission_->ECont()()();
487 
488  // Only include continuous phase absorption
490  absorptionEmission_->aCont()()();
491 
492  return a*G - E;
493 }
494 
495 
496 void Foam::radiation::fvDOM::updateBlackBodyEmission()
497 {
498  for (label j=0; j < nLambda_; j++)
499  {
500  blackBody_.correct(j, absorptionEmission_->bands(j));
501  }
502 }
503 
504 
506 {
507  G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
508  qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
509  qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
510  qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
511 
512  forAll(IRay_, rayI)
513  {
514  IRay_[rayI].addIntensity();
515  G_ += IRay_[rayI].I()*IRay_[rayI].omega();
516  qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
517  qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
518  qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
519  }
520 }
521 
522 
524 (
525  const word& name,
526  label& rayId,
527  label& lambdaId
528 ) const
529 {
530  // assuming name is in the form: CHARS_rayId_lambdaId
531  size_type i1 = name.find_first_of("_");
532  size_type i2 = name.find_last_of("_");
533 
534  rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
535  lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
536 }
537 
538 
539 // ************************************************************************* //
Collection of constants.
const volScalarField & a() const
Const access to total absorption coefficient.
Definition: fvDOM.H:69
dictionary coeffs_
Radiation model dictionary.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:431
#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].
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:644
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
Macros for easy insertion into run-time selection tables.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
void setRayIdLambdaId(const word &name, label &rayId, label &lambdaId) const
Set the rayId and lambdaId from by decomposing an intensity.
Definition: fvDOM.C:524
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:505
const fvMesh & mesh_
Reference to the mesh database.
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
Top level model for radiation modelling.
word timeName
Definition: getTimeIndex.H:3
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
label readLabel(Istream &is)
Definition: label.H:64
const volScalarField & G() const
Const access to incident radiation field.
Definition: fvDOM.H:84
bool read()
Read radiation properties dictionary.
Definition: fvDOM.C:397
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
dimensionedScalar pow3(const dimensionedScalar &ds)
virtual tmp< volScalarField::Internal > Ru() const
Source term component (constant)
Definition: fvDOM.C:479
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:170
virtual ~fvDOM()
Destructor.
Definition: fvDOM.C:391
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:419
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:53
const scalar piByTwo(0.5 *pi)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
virtual tmp< volScalarField > Rp() const
Source term component (for power of T^4)
Definition: fvDOM.C:456
#define addToRadiationRunTimeSelectionTables(model)
autoPtr< absorptionEmissionModel > absorptionEmission_
Absorption/emission model.
bool found
autoPtr< radiation::radiationModel > radiation(radiation::radiationModel::New(T))
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:576