radiativeIntensityRay.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 "radiativeIntensityRay.H"
27 #include "fvm.H"
28 #include "fvDOM.H"
29 #include "constants.H"
30 
31 using namespace Foam::constant;
32 
34 (
35  "ILambda"
36 );
37 
38 
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
42 (
43  const fvDOM& dom,
44  const fvMesh& mesh,
45  const scalar phi,
46  const scalar theta,
47  const scalar deltaPhi,
48  const scalar deltaTheta,
49  const label nLambda,
50  const absorptionEmissionModel& absorptionEmission,
51  const blackBodyEmission& blackBody,
52  const label rayId
53 )
54 :
55  dom_(dom),
56  mesh_(mesh),
57  absorptionEmission_(absorptionEmission),
58  blackBody_(blackBody),
59  I_
60  (
61  IOobject
62  (
63  "I" + name(rayId),
64  mesh_.time().timeName(),
65  mesh_,
68  ),
69  mesh_,
71  ),
72  qr_
73  (
74  IOobject
75  (
76  "qr" + name(rayId),
77  mesh_.time().timeName(),
78  mesh_,
81  ),
82  mesh_,
84  ),
85  qin_
86  (
87  IOobject
88  (
89  "qin" + name(rayId),
90  mesh_.time().timeName(),
91  mesh_,
94  ),
95  mesh_,
97  ),
98  qem_
99  (
100  IOobject
101  (
102  "qem" + name(rayId),
103  mesh_.time().timeName(),
104  mesh_,
107  ),
108  mesh_,
110  ),
111  d_(Zero),
112  dAve_(Zero),
113  theta_(theta),
114  phi_(phi),
115  omega_(0.0),
116  nLambda_(nLambda),
117  ILambda_(nLambda),
118  myRayId_(rayId)
119 {
120  scalar sinTheta = Foam::sin(theta);
121  scalar cosTheta = Foam::cos(theta);
122  scalar sinPhi = Foam::sin(phi);
123  scalar cosPhi = Foam::cos(phi);
124 
125  omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
126  d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
127  dAve_ = vector
128  (
129  sinPhi
130  *Foam::sin(0.5*deltaPhi)
131  *(deltaTheta - Foam::cos(2.0*theta)
132  *Foam::sin(deltaTheta)),
133  cosPhi
134  *Foam::sin(0.5*deltaPhi)
135  *(deltaTheta - Foam::cos(2.0*theta)
136  *Foam::sin(deltaTheta)),
137  0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
138  );
139 
140  // Transform directions so that they fall inside the bounds of reduced
141  // dimension cases
142  if (mesh_.nSolutionD() == 2)
143  {
144  vector meshDir(vector::zero);
145  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
146  {
147  if (mesh_.geometricD()[cmpt] == -1)
148  {
149  meshDir[cmpt] = 1;
150  }
151  }
152  const vector normal(vector(0, 0, 1));
153 
154  const tensor coordRot = rotationTensor(normal, meshDir);
155 
156  dAve_ = coordRot & dAve_;
157  d_ = coordRot & d_;
158  }
159  else if (mesh_.nSolutionD() == 1)
160  {
161  vector meshDir(vector::zero);
162  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
163  {
164  if (mesh_.geometricD()[cmpt] == 1)
165  {
166  meshDir[cmpt] = 1;
167  }
168  }
169  const vector normal(vector(1, 0, 0));
170 
171  dAve_ = (dAve_ & normal)*meshDir;
172  d_ = (d_ & normal)*meshDir;
173  }
174 
175 
176  autoPtr<volScalarField> IDefaultPtr;
177 
178  forAll(ILambda_, lambdaI)
179  {
180  IOobject IHeader
181  (
182  intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
183  mesh_.time().timeName(),
184  mesh_,
187  );
188 
189  // Check if field exists and can be read
190  if (IHeader.typeHeaderOk<volScalarField>(true))
191  {
192  ILambda_.set
193  (
194  lambdaI,
195  new volScalarField(IHeader, mesh_)
196  );
197  }
198  else
199  {
200  // Demand driven load the IDefault field
201  if (!IDefaultPtr.valid())
202  {
203  IDefaultPtr.reset
204  (
205  new volScalarField
206  (
207  IOobject
208  (
209  "IDefault",
210  mesh_.time().timeName(),
211  mesh_,
214  ),
215  mesh_
216  )
217  );
218  }
219 
220  // Reset the MUST_READ flag
221  IOobject noReadHeader(IHeader);
222  noReadHeader.readOpt() = IOobject::NO_READ;
223 
224  ILambda_.set
225  (
226  lambdaI,
227  new volScalarField(noReadHeader, IDefaultPtr())
228  );
229  }
230  }
231 }
232 
233 
234 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
235 
237 {}
238 
239 
240 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
241 
243 {
244  // Reset boundary heat flux to zero
245  qr_.boundaryFieldRef() = 0.0;
246 
247  scalar maxResidual = -great;
248 
249  const surfaceScalarField Ji(dAve_ & mesh_.Sf());
250 
251  forAll(ILambda_, lambdaI)
252  {
253  const volScalarField& k = dom_.aLambda(lambdaI);
254 
255  fvScalarMatrix IiEq
256  (
257  fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
258  + fvm::Sp(k*omega_, ILambda_[lambdaI])
259  ==
260  1.0/constant::mathematical::pi*omega_
261  *(
262  // Remove aDisp from k
263  (k - absorptionEmission_.aDisp(lambdaI))
264  *blackBody_.bLambda(lambdaI)
265 
266  + absorptionEmission_.E(lambdaI)/4
267  )
268  );
269 
270  IiEq.relax();
271 
272  const solverPerformance ILambdaSol = solve(IiEq, "Ii");
273 
274  const scalar initialRes =
275  ILambdaSol.initialResidual()*omega_/dom_.omegaMax();
276 
277  maxResidual = max(initialRes, maxResidual);
278  }
279 
280  return maxResidual;
281 }
282 
283 
285 {
287 
288  forAll(ILambda_, lambdaI)
289  {
290  I_ += ILambda_[lambdaI];
291  }
292 }
293 
294 
295 // ************************************************************************* //
Collection of constants.
Finite Volume Discrete Ordinates Method. Solves the RTE equation for n directions in a participating ...
Definition: fvDOM.H:74
#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 addIntensity()
Add radiative intensities from all the bands.
void reset(T *=nullptr)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
uint8_t direction
Definition: direction.H:45
scalar correct()
Update radiative intensity on i direction.
radiativeIntensityRay(const fvDOM &dom, const fvMesh &mesh, const scalar phi, const scalar theta, const scalar deltaPhi, const scalar deltaTheta, const label lambda, const absorptionEmissionModel &absEmmModel_, const blackBodyEmission &blackBody, const label rayId)
Construct form components.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
label k
Boltzmann constant.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from unit vector n1 to n2.
Definition: transform.H:47
T cosPhi(const Vector< T > &a, const Vector< T > &b, const T &tolerance=small)
Calculate angle between a and b in radians.
Definition: vectorTools.H:105
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
const dimensionSet dimTime
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:99
dimensionedScalar cos(const dimensionedScalar &ds)
const Type & initialResidual() const
Return initial residual.
A class for handling words, derived from string.
Definition: word.H:59
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics...
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
static const zero Zero
Definition: zero.H:97
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
dimensionedScalar sin(const dimensionedScalar &ds)
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:521
const dimensionSet dimMass
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
Model to supply absorption and emission coefficients for radiation modelling.
readOption readOpt() const
Definition: IOobject.H:353
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92