FreeFOAM The Cross-Platform CFD Toolkit
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) 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 "radiativeIntensityRay.H"
27 #include <finiteVolume/fvm.H>
28 #include <radiation/fvDOM.H>
29 
30 const Foam::word
32 
33 
34 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
35 
36 Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
37 (
38  const fvDOM& dom,
39  const fvMesh& mesh,
40  const scalar phi,
41  const scalar theta,
42  const scalar deltaPhi,
43  const scalar deltaTheta,
44  const label nLambda,
45  const absorptionEmissionModel& absorptionEmission,
46  const blackBodyEmission& blackBody,
47  const label rayId
48 )
49 :
50  dom_(dom),
51  mesh_(mesh),
52  absorptionEmission_(absorptionEmission),
53  blackBody_(blackBody),
54  I_
55  (
56  IOobject
57  (
58  "I" + name(rayId),
59  mesh_.time().timeName(),
60  mesh_,
61  IOobject::NO_READ,
62  IOobject::NO_WRITE
63  ),
64  mesh_,
66  ),
67  Qr_
68  (
69  IOobject
70  (
71  "Qr" + name(rayId),
72  mesh_.time().timeName(),
73  mesh_,
74  IOobject::NO_READ,
75  IOobject::NO_WRITE
76  ),
77  mesh_,
79  ),
80  Qin_
81  (
82  IOobject
83  (
84  "Qin" + name(rayId),
85  mesh_.time().timeName(),
86  mesh_,
87  IOobject::NO_READ,
88  IOobject::NO_WRITE
89  ),
90  mesh_,
92  ),
93  Qem_
94  (
95  IOobject
96  (
97  "Qem" + name(rayId),
98  mesh_.time().timeName(),
99  mesh_,
100  IOobject::NO_READ,
101  IOobject::NO_WRITE
102  ),
103  mesh_,
104  dimensionedScalar("Qem", dimMass/pow3(dimTime), 0.0)
105  ),
106  d_(vector::zero),
107  dAve_(vector::zero),
108  theta_(theta),
109  phi_(phi),
110  omega_(0.0),
111  nLambda_(nLambda),
112  ILambda_(nLambda)
113 {
114  scalar sinTheta = Foam::sin(theta);
115  scalar cosTheta = Foam::cos(theta);
116  scalar sinPhi = Foam::sin(phi);
117  scalar cosPhi = Foam::cos(phi);
118 
119  omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
120  d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
121  dAve_ = vector
122  (
123  sinPhi
124  *Foam::sin(0.5*deltaPhi)
125  *(deltaTheta - Foam::cos(2.0*theta)
126  *Foam::sin(deltaTheta)),
127  cosPhi
128  *Foam::sin(0.5*deltaPhi)
129  *(deltaTheta - Foam::cos(2.0*theta)
130  *Foam::sin(deltaTheta)),
131  0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
132  );
133 
134 
135  autoPtr<volScalarField> IDefaultPtr;
136 
137  forAll(ILambda_, lambdaI)
138  {
139  IOobject IHeader
140  (
141  intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
142  mesh_.time().timeName(),
143  mesh_,
144  IOobject::MUST_READ,
145  IOobject::AUTO_WRITE
146  );
147 
148  // check if field exists and can be read
149  if (IHeader.headerOk())
150  {
151  ILambda_.set
152  (
153  lambdaI,
154  new volScalarField(IHeader, mesh_)
155  );
156  }
157  else
158  {
159  // Demand driven load the IDefault field
160  if (!IDefaultPtr.valid())
161  {
162  IDefaultPtr.reset
163  (
164  new volScalarField
165  (
166  IOobject
167  (
168  "IDefault",
169  mesh_.time().timeName(),
170  mesh_,
171  IOobject::MUST_READ,
172  IOobject::NO_WRITE
173  ),
174  mesh_
175  )
176  );
177  }
178 
179  // Reset the MUST_READ flag
180  IOobject noReadHeader(IHeader);
181  noReadHeader.readOpt() = IOobject::NO_READ;
182 
183  ILambda_.set
184  (
185  lambdaI,
186  new volScalarField(noReadHeader, IDefaultPtr())
187  );
188  }
189  }
190 }
191 
192 
193 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
194 
196 {}
197 
198 
199 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
200 
202 {
203  // reset boundary heat flux to zero
204  Qr_.boundaryField() = 0.0;
205 
206  scalar maxResidual = -GREAT;
207 
208  forAll(ILambda_, lambdaI)
209  {
210  const volScalarField& k = dom_.aLambda(lambdaI);
211 
212  surfaceScalarField Ji = dAve_ & mesh_.Sf();
213 
214  fvScalarMatrix IiEq
215  (
216  fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
217  + fvm::Sp(k*omega_, ILambda_[lambdaI])
218  ==
220  *(
221  k*blackBody_.bLambda(lambdaI)
222  + absorptionEmission_.ECont(lambdaI)/4
223  )
224  );
225 
226  IiEq.relax();
227 
228  scalar eqnResidual = solve
229  (
230  IiEq,
231  mesh_.solver("Ii")
232  ).initialResidual();
233 
234  maxResidual = max(eqnResidual, maxResidual);
235  }
236 
237  return maxResidual;
238 }
239 
240 
242 {
243  I_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
244 
245  forAll(ILambda_, lambdaI)
246  {
247  I_ += absorptionEmission_.addIntensity(lambdaI, ILambda_[lambdaI]);
248  }
249 }
250 
251 
252 // ************************ vim: set sw=4 sts=4 et: ************************ //