FreeFOAM The Cross-Platform CFD Toolkit
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) 2008-2010 OpenCFD Ltd.
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 
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  namespace radiation
37  {
38  defineTypeNameAndDebug(fvDOM, 0);
39 
41  (
42  radiationModel,
43  fvDOM,
44  dictionary
45  );
46  }
47 }
48 
49 
50 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
51 
52 Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
53 :
54  radiationModel(typeName, T),
55  G_
56  (
57  IOobject
58  (
59  "G",
60  mesh_.time().timeName(),
61  mesh_,
62  IOobject::NO_READ,
63  IOobject::AUTO_WRITE
64  ),
65  mesh_,
67  ),
68  Qr_
69  (
70  IOobject
71  (
72  "Qr",
73  mesh_.time().timeName(),
74  mesh_,
75  IOobject::NO_READ,
76  IOobject::AUTO_WRITE
77  ),
78  mesh_,
80  ),
81  Qem_
82  (
83  IOobject
84  (
85  "Qem",
86  mesh_.time().timeName(),
87  mesh_,
88  IOobject::NO_READ,
89  IOobject::NO_WRITE
90  ),
91  mesh_,
93  ),
94  Qin_
95  (
96  IOobject
97  (
98  "Qin",
99  mesh_.time().timeName(),
100  mesh_,
101  IOobject::NO_READ,
102  IOobject::NO_WRITE
103  ),
104  mesh_,
105  dimensionedScalar("Qin", dimMass/pow3(dimTime), 0.0)
106  ),
107  a_
108  (
109  IOobject
110  (
111  "a",
112  mesh_.time().timeName(),
113  mesh_,
114  IOobject::NO_READ,
115  IOobject::NO_WRITE
116  ),
117  mesh_,
119  ),
120  e_
121  (
122  IOobject
123  (
124  "e",
125  mesh_.time().timeName(),
126  mesh_,
127  IOobject::NO_READ,
128  IOobject::NO_WRITE
129  ),
130  mesh_,
132  ),
133  E_
134  (
135  IOobject
136  (
137  "E",
138  mesh_.time().timeName(),
139  mesh_,
140  IOobject::NO_READ,
141  IOobject::NO_WRITE
142  ),
143  mesh_,
145  ),
146  nTheta_(readLabel(coeffs_.lookup("nTheta"))),
147  nPhi_(readLabel(coeffs_.lookup("nPhi"))),
148  nRay_(0),
149  nLambda_(absorptionEmission_->nBands()),
150  aLambda_(nLambda_),
151  blackBody_(nLambda_, T),
152  IRay_(0),
153  convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
154  maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50))
155 {
156  if (mesh_.nSolutionD() == 3) //3D
157  {
158  nRay_ = 4*nPhi_*nTheta_;
159  IRay_.setSize(nRay_);
160  scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
161  scalar deltaTheta = mathematicalConstant::pi/nTheta_;
162  label i = 0;
163  for (label n = 1; n <= nTheta_; n++)
164  {
165  for (label m = 1; m <= 4*nPhi_; m++)
166  {
167  scalar thetai = (2.0*n - 1.0)*deltaTheta/2.0;
168  scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
169  IRay_.set
170  (
171  i,
173  (
174  *this,
175  mesh_,
176  phii,
177  thetai,
178  deltaPhi,
179  deltaTheta,
180  nLambda_,
182  blackBody_,
183  i
184  )
185  );
186  i++;
187  }
188  }
189  }
190  else
191  {
192  if (mesh_.nSolutionD() == 2) //2D (X & Y)
193  {
194  scalar thetai = mathematicalConstant::piByTwo;
195  scalar deltaTheta = mathematicalConstant::pi;
196  nRay_ = 4*nPhi_;
197  IRay_.setSize(nRay_);
198  scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
199  label i = 0;
200  for (label m = 1; m <= 4*nPhi_; m++)
201  {
202  scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
203  IRay_.set
204  (
205  i,
207  (
208  *this,
209  mesh_,
210  phii,
211  thetai,
212  deltaPhi,
213  deltaTheta,
214  nLambda_,
216  blackBody_,
217  i
218  )
219  );
220  i++;
221  }
222  }
223  else //1D (X)
224  {
225  scalar thetai = mathematicalConstant::piByTwo;
226  scalar deltaTheta = mathematicalConstant::pi;
227  nRay_ = 2;
228  IRay_.setSize(nRay_);
229  scalar deltaPhi = mathematicalConstant::pi;
230  label i = 0;
231  for (label m = 1; m <= 2; m++)
232  {
233  scalar phii = (2.0*m - 1.0)*deltaPhi/2.0;
234  IRay_.set
235  (
236  i,
238  (
239  *this,
240  mesh_,
241  phii,
242  thetai,
243  deltaPhi,
244  deltaTheta,
245  nLambda_,
247  blackBody_,
248  i
249  )
250  );
251  i++;
252  }
253 
254  }
255  }
256 
257 
258  // Construct absorption field for each wavelength
259  forAll(aLambda_, lambdaI)
260  {
261  aLambda_.set
262  (
263  lambdaI,
264  new volScalarField
265  (
266  IOobject
267  (
268  "aLambda_" + Foam::name(lambdaI) ,
269  mesh_.time().timeName(),
270  mesh_,
273  ),
274  a_
275  )
276  );
277  }
278 
279  Info<< "fvDOM : Allocated " << IRay_.size()
280  << " rays with average orientation:" << nl;
281  forAll (IRay_, i)
282  {
283  Info<< '\t' << IRay_[i].I().name()
284  << '\t' << IRay_[i].dAve() << nl;
285  }
286  Info<< endl;
287 }
288 
289 
290 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
291 
293 {}
294 
295 
296 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
297 
299 {
300  if (radiationModel::read())
301  {
302 // Only reading solution parameters - not changing ray geometry
303 
304  coeffs_.readIfPresent("convergence", convergence_);
305  coeffs_.readIfPresent("maxIter", maxIter_);
306 
307  return true;
308  }
309  else
310  {
311  return false;
312  }
313 }
314 
315 
317 {
318  absorptionEmission_->correct(a_, aLambda_);
319 
320  updateBlackBodyEmission();
321 
322  scalar maxResidual = 0.0;
323  label radIter = 0;
324  do
325  {
326  radIter++;
327  forAll(IRay_, rayI)
328  {
329  maxResidual = 0.0;
330  scalar maxBandResidual = IRay_[rayI].correct();
331  maxResidual = max(maxBandResidual, maxResidual);
332  }
333 
334  Info << "Radiation solver iter: " << radIter << endl;
335 
336  } while(maxResidual > convergence_ && radIter < maxIter_);
337 
338  updateG();
339 }
340 
341 
343 {
344  return tmp<volScalarField>
345  (
346  new volScalarField
347  (
348  IOobject
349  (
350  "Rp",
351  mesh_.time().timeName(),
352  mesh_,
355  false
356  ),
357  4.0*a_*radiation::sigmaSB //absorptionEmission_->a()
358  )
359  );
360 }
361 
362 
365 {
366 
368  G_.dimensionedInternalField();
370  absorptionEmission_->ECont()().dimensionedInternalField();
372  a_.dimensionedInternalField(); //absorptionEmission_->aCont()()
373 
374  return a*G - E;
375 }
376 
377 
378 void Foam::radiation::fvDOM::updateBlackBodyEmission()
379 {
380  for (label j=0; j < nLambda_; j++)
381  {
382  blackBody_.correct(j, absorptionEmission_->bands(j));
383  }
384 }
385 
386 
388 {
389  G_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
390  Qr_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
391  Qem_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
392  Qin_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
393 
394  forAll(IRay_, rayI)
395  {
396  IRay_[rayI].addIntensity();
397  G_ += IRay_[rayI].I()*IRay_[rayI].omega();
398  Qr_.boundaryField() += IRay_[rayI].Qr().boundaryField();
399  Qem_.boundaryField() += IRay_[rayI].Qem().boundaryField();
400  Qin_.boundaryField() += IRay_[rayI].Qin().boundaryField();
401  }
402 }
403 
404 
406 (
407  const word& name,
408  label& rayId,
409  label& lambdaId
410 ) const
411 {
412  // assuming name is in the form: CHARS_rayId_lambdaId
413  size_type i1 = name.find_first_of("_");
414  size_type i2 = name.find_last_of("_");
415 
416  rayId = readLabel(IStringStream(name.substr(i1+1, i2-1))());
417  lambdaId = readLabel(IStringStream(name.substr(i2+1, name.size()-1))());
418 }
419 
420 
421 // ************************ vim: set sw=4 sts=4 et: ************************ //