radiativeIntensityRay.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 "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 
41 Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
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_,
109  dimensionedScalar("qem", dimMass/pow3(dimTime), 0.0)
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  autoPtr<volScalarField> IDefaultPtr;
141 
142  forAll(ILambda_, lambdaI)
143  {
144  IOobject IHeader
145  (
146  intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
147  mesh_.time().timeName(),
148  mesh_,
151  );
152 
153  // Check if field exists and can be read
154  if (IHeader.typeHeaderOk<volScalarField>(true))
155  {
156  ILambda_.set
157  (
158  lambdaI,
159  new volScalarField(IHeader, mesh_)
160  );
161  }
162  else
163  {
164  // Demand driven load the IDefault field
165  if (!IDefaultPtr.valid())
166  {
167  IDefaultPtr.reset
168  (
169  new volScalarField
170  (
171  IOobject
172  (
173  "IDefault",
174  mesh_.time().timeName(),
175  mesh_,
178  ),
179  mesh_
180  )
181  );
182  }
183 
184  // Reset the MUST_READ flag
185  IOobject noReadHeader(IHeader);
186  noReadHeader.readOpt() = IOobject::NO_READ;
187 
188  ILambda_.set
189  (
190  lambdaI,
191  new volScalarField(noReadHeader, IDefaultPtr())
192  );
193  }
194  }
195 }
196 
197 
198 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
199 
201 {}
202 
203 
204 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
205 
207 {
208  // Reset boundary heat flux to zero
209  qr_.boundaryFieldRef() = 0.0;
210 
211  scalar maxResidual = -GREAT;
212 
213  forAll(ILambda_, lambdaI)
214  {
215  const volScalarField& k = dom_.aLambda(lambdaI);
216 
217  const surfaceScalarField Ji(dAve_ & mesh_.Sf());
218 
219  fvScalarMatrix IiEq
220  (
221  fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
222  + fvm::Sp(k*omega_, ILambda_[lambdaI])
223  ==
224  1.0/constant::mathematical::pi*omega_
225  *(
226  // Remove aDisp from k
227  (k - absorptionEmission_.aDisp(lambdaI))
228  *blackBody_.bLambda(lambdaI)
229 
230  + absorptionEmission_.E(lambdaI)/4
231  )
232  );
233 
234  IiEq.relax();
235 
236  const solverPerformance ILambdaSol = solve
237  (
238  IiEq,
239  mesh_.solver("Ii")
240  );
241 
242  const scalar initialRes =
243  ILambdaSol.initialResidual()*omega_/dom_.omegaMax();
244 
245  maxResidual = max(initialRes, maxResidual);
246  }
247 
248  return maxResidual;
249 }
250 
251 
253 {
254  I_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
255 
256  forAll(ILambda_, lambdaI)
257  {
258  I_ += ILambda_[lambdaI];
259  }
260 }
261 
262 
263 // ************************************************************************* //
Collection of constants.
Class black body emission.
#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
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
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 > &)
Model to supply absorption and emission coefficients for radiation modelling.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
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
scalar correct()
Update radiative intensity on i direction.
void reset(T *=0)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
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< Type > solve(fvMatrix< Type > &, const dictionary &)
Solve returning the solution statistics given convergence tolerance.
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:91
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 addIntensity()
Add radiative intensities from all the bands.
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
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
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
Finite Volume Discrete Ordinates Method. Solves the RTE equation for n directions in a participating ...
Definition: fvDOM.H:80