radiativeIntensityRay.C 8.21 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
6
     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd
7
8
9
10
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11
12
13
14
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
15
16
17
18
19
20
21

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    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
27
28

\*---------------------------------------------------------------------------*/

#include "radiativeIntensityRay.H"
#include "fvm.H"
#include "fvDOM.H"
mattijs's avatar
mattijs committed
29
30
31
#include "constants.H"

using namespace Foam::constant;
32

33
34
const Foam::word
Foam::radiation::radiativeIntensityRay::intensityPrefix("ILambda");
35
36


37
38
39
40
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
(
41
42
    const fvDOM& dom,
    const fvMesh& mesh,
43
44
45
46
47
    const scalar phi,
    const scalar theta,
    const scalar deltaPhi,
    const scalar deltaTheta,
    const label nLambda,
Andrew Heather's avatar
tidying    
Andrew Heather committed
48
    const absorptionEmissionModel& absorptionEmission,
49
50
    const blackBodyEmission& blackBody,
    const label rayId
51
52
)
:
53
    dom_(dom),
54
    mesh_(mesh),
Andrew Heather's avatar
tidying    
Andrew Heather committed
55
    absorptionEmission_(absorptionEmission),
56
57
58
59
60
61
62
63
64
65
66
67
    blackBody_(blackBody),
    I_
    (
        IOobject
        (
            "I" + name(rayId),
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
68
        dimensionedScalar(dimMass/pow3(dimTime), Zero)
69
    ),
70
    qr_
71
72
73
    (
        IOobject
        (
74
            "qr" + name(rayId),
75
76
77
78
79
80
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
81
        dimensionedScalar(dimMass/pow3(dimTime), Zero)
82
    ),
83
    qin_
84
85
86
    (
        IOobject
        (
87
            "qin" + name(rayId),
88
89
90
91
92
93
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
94
        dimensionedScalar(dimMass/pow3(dimTime), Zero)
95
    ),
96
    qem_
97
98
99
    (
        IOobject
        (
100
            "qem" + name(rayId),
101
102
103
104
105
106
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
107
        dimensionedScalar(dimMass/pow3(dimTime), Zero)
108
    ),
Henry Weller's avatar
Henry Weller committed
109
110
    d_(Zero),
    dAve_(Zero),
111
112
113
114
    theta_(theta),
    phi_(phi),
    omega_(0.0),
    nLambda_(nLambda),
Sergio Ferraris's avatar
Sergio Ferraris committed
115
116
    ILambda_(nLambda),
    myRayId_(rayId)
117
{
118
119
120
121
122
123
124
125
    scalar sinTheta = Foam::sin(theta);
    scalar cosTheta = Foam::cos(theta);
    scalar sinPhi = Foam::sin(phi);
    scalar cosPhi = Foam::cos(phi);

    omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
    d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
    dAve_ = vector
126
    (
127
128
        sinPhi
       *Foam::sin(0.5*deltaPhi)
129
       *(deltaTheta - Foam::cos(2.0*theta)
130
131
132
       *Foam::sin(deltaTheta)),
        cosPhi
       *Foam::sin(0.5*deltaPhi)
133
       *(deltaTheta - Foam::cos(2.0*theta)
134
       *Foam::sin(deltaTheta)),
135
        0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
136
137
    );

138
139
    if (mesh_.nSolutionD() == 2)
    {
140
141
142
143
144
145
146
147
148
149
150
        // Omega for 2D
        omega_ = deltaPhi;

        // dAve for 2D
        dAve_ = vector
        (
            2*sinPhi*Foam::sin(0.5*deltaPhi),
            2*cosPhi*Foam::sin(0.5*deltaPhi),
            0
        );

Andrew Heather's avatar
Andrew Heather committed
151
        vector meshDir(Zero);
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
        if (dom_.meshOrientation() != vector::zero)
        {
            meshDir = dom_.meshOrientation();
        }
        else
        {
            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
            {
                if (mesh_.geometricD()[cmpt] == -1)
                {
                    meshDir[cmpt] = 1;
                }
            }
        }
        const vector normal(vector(0, 0, 1));

        const tensor coordRot = rotationTensor(normal, meshDir);

        dAve_ = coordRot & dAve_;
        d_ = coordRot & d_;
172

173
174
175
    }
    else if (mesh_.nSolutionD() == 1)
    {
Andrew Heather's avatar
Andrew Heather committed
176
        vector meshDir(Zero);
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
        if (dom_.meshOrientation() != vector::zero)
        {
            meshDir = dom_.meshOrientation();
        }
        else
        {
            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
            {
                if (mesh_.geometricD()[cmpt] == 1)
                {
                    meshDir[cmpt] = 1;
                }
            }
        }
        const vector normal(vector(1, 0, 0));

        dAve_ = (dAve_ & normal)*meshDir;
        d_ = (d_ & normal)*meshDir;
195
196
197

        // Omega normalization for 1D
        omega_ /= 2;
198
199
    }

mattijs's avatar
mattijs committed
200
201
    autoPtr<volScalarField> IDefaultPtr;

Andrew Heather's avatar
Andrew Heather committed
202
    forAll(ILambda_, lambdaI)
203
    {
204
        IOobject IHeader
205
        (
206
            intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
207
208
            mesh_.time().timeName(),
            mesh_,
209
210
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
211
        );
212

Andrew Heather's avatar
Andrew Heather committed
213
        // Check if field exists and can be read
214
        if (IHeader.typeHeaderOk<volScalarField>(true))
215
        {
Andrew Heather's avatar
Andrew Heather committed
216
            ILambda_.set
217
            (
Andrew Heather's avatar
Andrew Heather committed
218
                lambdaI,
219
                new volScalarField(IHeader, mesh_)
220
221
222
223
            );
        }
        else
        {
mattijs's avatar
mattijs committed
224
225
226
227
            // Demand driven load the IDefault field
            if (!IDefaultPtr.valid())
            {
                IDefaultPtr.reset
228
                (
mattijs's avatar
mattijs committed
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
                    new volScalarField
                    (
                        IOobject
                        (
                            "IDefault",
                            mesh_.time().timeName(),
                            mesh_,
                            IOobject::MUST_READ,
                            IOobject::NO_WRITE
                        ),
                        mesh_
                    )
                );
            }

            // Reset the MUST_READ flag
            IOobject noReadHeader(IHeader);
            noReadHeader.readOpt() = IOobject::NO_READ;
247

Andrew Heather's avatar
Andrew Heather committed
248
            ILambda_.set
249
            (
Andrew Heather's avatar
Andrew Heather committed
250
                lambdaI,
mattijs's avatar
mattijs committed
251
                new volScalarField(noReadHeader, IDefaultPtr())
252
253
254
255
256
            );
        }
    }
}

257

258
259
260
261
262
263
264
265
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::radiation::radiativeIntensityRay::~radiativeIntensityRay()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

266
Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
267
{
268
    // Reset boundary heat flux to zero
269
    qr_.boundaryFieldRef() = 0.0;
270

271
    scalar maxResidual = -GREAT;
272

Andrew Heather's avatar
Andrew Heather committed
273
    forAll(ILambda_, lambdaI)
274
    {
Andrew Heather's avatar
Andrew Heather committed
275
        const volScalarField& k = dom_.aLambda(lambdaI);
276

277
        const surfaceScalarField Ji(dAve_ & mesh_.Sf());
278

279
280
281
282
283
284
285
286
287
288
289
290
291
        fvScalarMatrix IiEq
        (
            fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
          + fvm::Sp(k*omega_, ILambda_[lambdaI])
        ==
            1.0/constant::mathematical::pi*omega_
           *(
                (k - absorptionEmission_.aDisp(lambdaI))
               *blackBody_.bLambda(lambdaI)

              + absorptionEmission_.E(lambdaI)/4
            )
        );
292

293
        IiEq.relax();
Sergio Ferraris's avatar
Sergio Ferraris committed
294
295

        const solverPerformance ILambdaSol = solve
296
        (
297
            IiEq,
298
            mesh_.solver("Ii")
Sergio Ferraris's avatar
Sergio Ferraris committed
299
300
301
302
        );

        const scalar initialRes =
            ILambdaSol.initialResidual()*omega_/dom_.omegaMax();
303

Sergio Ferraris's avatar
Sergio Ferraris committed
304
        maxResidual = max(initialRes, maxResidual);
305
    }
306

307
308
309
    return maxResidual;
}

310

311
312
void Foam::radiation::radiativeIntensityRay::addIntensity()
{
313
    I_ = dimensionedScalar(dimMass/pow3(dimTime), Zero);
314

Andrew Heather's avatar
Andrew Heather committed
315
    forAll(ILambda_, lambdaI)
316
    {
317
        I_ += ILambda_[lambdaI];
318
319
320
    }
}

321

322
// ************************************************************************* //