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-2022 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().name(),
65  mesh_,
66  IOobject::NO_READ,
67  IOobject::NO_WRITE
68  ),
69  mesh_,
71  ),
72  qr_
73  (
74  IOobject
75  (
76  "qr" + name(rayId),
77  mesh_.time().name(),
78  mesh_,
79  IOobject::NO_READ,
80  IOobject::NO_WRITE
81  ),
82  mesh_,
84  ),
85  qin_
86  (
87  IOobject
88  (
89  "qin" + name(rayId),
90  mesh_.time().name(),
91  mesh_,
92  IOobject::NO_READ,
93  IOobject::NO_WRITE
94  ),
95  mesh_,
97  ),
98  qem_
99  (
100  IOobject
101  (
102  "qem" + name(rayId),
103  mesh_.time().name(),
104  mesh_,
105  IOobject::NO_READ,
106  IOobject::NO_WRITE
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  {
181  (
182  intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
183  mesh_.time().name(),
184  mesh_,
187  );
188 
189  // Check if field exists and can be read
190  if (IHeader.headerOk())
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().name(),
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 // ************************************************************************* //
label k
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Generic GeometricField class.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
readOption readOpt() const
Definition: IOobject.H:360
SolverPerformance is the class returned by the LduMatrix solver containing performance statistics.
const Type & initialResidual() const
Return initial residual.
static const Form zero
Definition: VectorSpace.H:118
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:105
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
void reset(T *=nullptr)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
const word & name() const
Return const reference to name.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:603
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:418
label nSolutionD() const
Return the number of valid solved-for dimensions in the mesh.
Definition: polyMesh.C:1017
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:989
Model to supply absorption and emission coefficients for radiation modelling.
Finite Volume Discrete Ordinates Method. Solves the RTE equation for n directions in a participating ...
Definition: fvDOM.H:77
scalar theta() const
Return the theta angle.
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.
scalar phi() const
Return the phi angle.
void addIntensity()
Add radiative intensities from all the bands.
Templated form of IOobject providing type information for file reading and header type checking.
Definition: IOobject.H:531
bool headerOk()
Read header (uses typeGlobalFile to find file) and check.
A class for handling words, derived from string.
Definition: word.H:62
Collection of constants.
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const VolField< Type > &)
static const zero Zero
Definition: zero.H:97
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from unit vector n1 to n2.
Definition: transform.H:47
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
dimensionedScalar pow3(const dimensionedScalar &ds)
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionSet dimTime
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:64
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const dimensionSet dimMass
dimensionedScalar cos(const dimensionedScalar &ds)
uint8_t direction
Definition: direction.H:45
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.