fvDOM Model

To understand the how fvDOM model is used and implemented in OpenFOAM, we use tutorial smallPoolFire3D as an example. We have discussed the p1, which is also a radiation model implemented in OpenFOAM. Here, we will discuss another radiation model used in OpenFOAM called fvDOM (Finite Volume Discrete Ordinates Method). The implementation logic of this model in OpenFOAM is quite similar to how we did in P1 model. Here we will not repeat what we have explained in p1 model.

Before starting, the the mathematical form of the fvDOM model in ANSYS FLUENT will be shown:

In OpenFOAM, the last term on the rhs of Eq.() is not considered. The fvDOM equation in OpenFOAM will be shown later. Now, let’s go through the code of tutorial smallPoolFire3D to see how fireFoam uses the fvDOM model. Everything works as same as coalCombustionFoam until we choose the model. So, here I will directly go to the model selection part. In Foam::radiationModel::New, we will search for the models we input in the radiationProperties files:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
radiation on;

radiationModel  fvDOM;

fvDOMCoeffs
{
    nPhi        2;      // azimuthal angles in PI/2 on X-Y.(from Y to X)
    nTheta      2;      // polar angles in PI (from Z to X-Y plane)
    tolerance   1e-1;   // convergence tolerance for radiation iteration
    maxIter     1;      // maximum number of iterations
}

// Number of flow iterations per radiation iteration
solverFreq 10;

absorptionEmissionModel greyMeanCombustion;

constantCoeffs
{
    absorptivity    0.01;
    emissivity      0.01;
    E               0;
}

greyMeanCombustionCoeffs
{
    lookUpTableFileName      none;

    EhrrCoeff                0.0;

    CO2
    {
        Tcommon         200;   // Common Temp
        invTemp         true;   // Is the polynomio using inverse temperature.
        Tlow            200;   // Low Temp
        Thigh           2500;  // High Temp

        loTcoeffs       // coefss for T < Tcommon
        (
            0           //  a0            +
            0           //  a1*T          +
            0           //  a2*T^(+/-)2   +
            0           //  a3*T^(+/-)3   +
            0           //  a4*T^(+/-)4   +
            0           //  a5*T^(+/-)5   +
        );
        hiTcoeffs        // coefss for T > Tcommon
        (
            18.741
            -121.31e3
            273.5e6
            -194.05e9
            56.31e12
            -5.8169e15
        );

    }

    // ...

scatterModel    none;

sootModel       none;

From radiationModel fvDOM;, we know that fvDOM model is used. Then we will go to the constructor of class fvDOM in :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 Foam::radiationModels::fvDOM::fvDOM(const volScalarField& T)
 :
     radiationModel(typeName, T),
// ommit the initialization of G_, qr_, qem_, qin_, a_ here
     nTheta_(readLabel(coeffs_.lookup("nTheta"))),
     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
     nRay_(0),
     nLambda_(absorptionEmission_->nBands()),
     aLambda_(nLambda_),
     blackBody_(nLambda_, T),
     IRay_(0),
     tolerance_
     (
         coeffs_.found("convergence")
       ? readScalar(coeffs_.lookup("convergence"))
       : coeffs_.lookupOrDefault<scalar>("tolerance", 0)
     ),
     maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
     omegaMax_(0)
 {
     initialise();
 }

In the initialization list, the radiationModel constructor will be called first. typeName will be fvDOM here:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Foam::radiationModel::radiationModel(const word& type, const volScalarField& T)
:
    IOdictionary(createIOobject(T.mesh())),
    mesh_(T.mesh()),
    time_(T.time()),
    T_(T),
    coeffs_(subOrEmptyDict(type + "Coeffs")),
    solverFreq_(1),
    firstIter_(true),
    absorptionEmission_(nullptr),
    scatter_(nullptr),
    soot_(nullptr)
{
    initialise();
}

The coeffs_ here will be:

1
2
3
4
5
6
{
    nPhi            2;
    nTheta          2;
    tolerance       0.1;
    maxIter         1;
}

Then, we call initialise() to initialize the absorptionEmission_, scatter_ and soot_. This has already been explained in P1 and will not be explained again here. After calling the constructor of the radiationModel class, other members in the fvDOM class will be initialized with a zero value first. nTheta_ which is and nPhi_ which is , are number angles in and direction. They are user defined value as shown in coeffs_ above. In the figure below you will see how and are defined in the coordinate system. and will determine how many radiation rays we have and how many radiation transport equations (RTE) we are going to solve.

Angular Coordinate System

Let’s go back to the fvDOM constructor, nRay_ means the number of rays we are going to solve. This is determined by how we discretized the RTE using the DOM (discrete-ordinate method). Because theoretically, there should be infinite amount of rays. It means we have to solve infinite number of RTEs, which is indeed impossible. Therefore the situation is simplified to finite number of directions of rays with the same amount of RTEs. nLambda_ is the number of wavelength bands. aLambda_ is the wavelength total absorption coefficient, blackBody_ is the blackbody radiation intensity, IRay_ is the ray radiation intensity.

After the initialization in the fvDOM constructor, the function initialise() will be called. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
void Foam::radiationModels::fvDOM::initialise()
{
    // 3D
    if (mesh_.nSolutionD() == 3)
    {
        nRay_ = 4*nPhi_*nTheta_;
        IRay_.setSize(nRay_);
        const scalar deltaPhi = pi/(2.0*nPhi_);
        const scalar deltaTheta = pi/nTheta_;
        label i = 0;
        for (label n = 1; n <= nTheta_; n++)
        {
            for (label m = 1; m <= 4*nPhi_; m++)
            {
                const scalar thetai = (2*n - 1)*deltaTheta/2.0;
                const scalar phii = (2*m - 1)*deltaPhi/2.0;
                IRay_.set
                (
                    i,
                    new radiativeIntensityRay
                    (
                        *this,
                        mesh_,
                        phii,
                        thetai,
                        deltaPhi,
                        deltaTheta,
                        nLambda_,
                        absorptionEmission_,
                        blackBody_,
                        i
                    )
                );
                i++;
            }
        }
    }
    //2D
    // ...
    //1D
    //...
    
    // Construct absorption field for each wavelength
    forAll(aLambda_, lambdaI)
    {
        aLambda_.set
        (
            lambdaI,
            new volScalarField
            (
                IOobject
                (
                    "aLambda_" + Foam::name(lambdaI) ,
                    mesh_.time().timeName(),
                    mesh_,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                a_
            )
        );
    }


    // Calculate the maximum solid angle
    forAll(IRay_, rayId)
    {
        if (omegaMax_ <  IRay_[rayId].omega())
        {
            omegaMax_ = IRay_[rayId].omega();
        }
    }


    Info<< typeName << ": Created " << IRay_.size() << " rays with average "
        << "directions (dAve) and solid angles (omega)" << endl;
    Info<< incrIndent;
    forAll(IRay_, rayId)
    {
        Info<< indent
            << "Ray " << IRay_[rayId].I().name() << ": "
            << "dAve = " << IRay_[rayId].dAve() << ", "
            << "omega = " << IRay_[rayId].omega() << endl;
    }
    Info<< decrIndent << endl;
}

We will only discuss the general 3D scenario here. The number of rays nRay_, is determined number of angles we what to discretize defined by nPhi_ and nTheta_. nPhi_ and nTheta_ are defined in the radiationProperties file. IRay_ is a pointer list which save informations of radiation rays. Here, the size of this pointer list is initialized by the number of rays as setSize(nRay_). Then, the increment variables deltaPhi and deltaTheta are defined. Here as the spherical coordinate system is used, the increment of and should also follow the rule of the spherical coordinate system. When we do the integral over the sphere, increases from 0 to and increases from 0 to . Here nPhi is azimuthal angles in on X-Y, nTheta is polar angles in (from Z to X-Y plane). Therefore, e.g. when nPhi and nTheta are both equal to one, there will be 1 and 4 . In this case, is equal to and will increase in a row like , , , and . Like what is described in the code:

1
2
3
4
5
6
7
8
9
10
for (label n = 1; n <= nTheta_; n++)
{
    for (label m = 1; m <= 4*nPhi_; m++)
    {
        const scalar thetai = (2*n - 1)*deltaTheta/2.0;
        const scalar phii = (2*m - 1)*deltaPhi/2.0;
        //.....
    }
    //.....
}

Then, inside the for loop, the IRay_ will be initialized using set(). Here i will count for the index of the radiation ray. new radiativeIntensityRay will call the constructor of the radiativeIntensityRay class. The radiativeIntensityRay constructor is very long, we will discuss it in a separate way. Firstly, the parameters the constructor takes and what we input in fvDOM.C:

1
2
3
4
5
6
7
8
9
10
11
12
13
Foam::radiationModels::radiativeIntensityRay::radiativeIntensityRay
(
    const fvDOM& dom,
    const fvMesh& mesh,
    const scalar phi,
    const scalar theta,
    const scalar deltaPhi,
    const scalar deltaTheta,
    const label nLambda,
    const absorptionEmissionModel& absorptionEmission,
    const blackBodyEmission& blackBody,
    const label rayId
)

and

1
2
3
4
5
6
7
8
9
10
11
12
13
new radiativeIntensityRay
(
    *this,
    mesh_,
    phii,
    thetai,
    deltaPhi,
    deltaTheta,
    nLambda_,
    absorptionEmission_,
    blackBody_,
    i
)

Then, dom_, mesh_, absorptionEmission_, and blackBody_ will be initialized with the input parameters. I_, qr_, qem_, and qem_ will be initialized with zero value. The radiation ray direction d_ and the average direction vector inside the solid angle dAve are initialized as zero. No more words about the constructor, please see the code below:

1
2
3
4
5
6
7
8
d_(Zero),
dAve_(Zero),
theta_(theta),
phi_(phi),
omega_(0.0),// the soild angle
nLambda_(nLambda), // number of wavelength
ILambda_(nLambda), //radiation intensity pointer list in certain wavelength
myRayId_(rayId)// index of the radiation ray

To save some time:

1
2
3
4
scalar sinTheta = Foam::sin(theta);
scalar cosTheta = Foam::cos(theta);
scalar sinPhi = Foam::sin(phi);
scalar cosPhi = Foam::cos(phi);

Then, the solid angle for “this” radiation ray will be calculated:

1
omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;

To explain this, some derivations will be shown here. In spherical coordinate system, the solid angle is defined as:

In one radiation ray direction, to calculate the solid angle, we need to do a integral over its angular range:

Then, the the direction of the ray in a unit sphere is defined:

1
d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);

In OpenFOAM, as we mentioned in radiationProperties, the angle is define from y->x but not x->y. That’s the reason why the d_ vector is defined as above. The average direction vector inside the solid angle is defined as:

1
2
3
4
5
6
7
8
9
10
11
12
dAve_ = vector
(
    sinPhi
   *Foam::sin(0.5*deltaPhi)
   *(deltaTheta - Foam::cos(2.0*theta)
   *Foam::sin(deltaTheta)),
    cosPhi
   *Foam::sin(0.5*deltaPhi)
   *(deltaTheta - Foam::cos(2.0*theta)
   *Foam::sin(deltaTheta)),
    0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
);

We will do the derivation here to explain how this vector is shown as this form. To solve the RTE using the fvDOM (discrete ordinate method), we need to do a finite volume integral over the spherical coordinate as what we did in Cartesian system. The RTE we are going to solve is:

As I mentioned before, OpenFOAM does not consider the last term on the rhs. So, the equation we are going to solve in OpenFOAM becomes:

It should be mentioned here, the is exactly the d_ (ray direction) as defined in the code. Now, the finite volume method will be applied to the RTE:

The last two terms on the rhs becomes:

and

The derivation of the first term on the lhs is a little bit different:

,where

For each term,

After the integration, the becomes:

This is defined as dAve_ in the radiativeIntensityRay constructor body:

1
2
3
4
5
6
7
8
9
10
11
12
dAve_ = vector
(
    sinPhi
   *Foam::sin(0.5*deltaPhi)
   *(deltaTheta - Foam::cos(2.0*theta)
   *Foam::sin(deltaTheta)),
    cosPhi
   *Foam::sin(0.5*deltaPhi)
   *(deltaTheta - Foam::cos(2.0*theta)
   *Foam::sin(deltaTheta)),
    0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
)

After the definition of dAve_, projections to 2D and 1D case are introduced. Here we will only discuss the 3D case. The initialization of ILambda_ will not be discussed here. When the for loop for different rays in the initialise() function is finished, the absorption coefficient for different wave length will be initialised:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// Construct absorption field for each wavelength
forAll(aLambda_, lambdaI)
{
    aLambda_.set
    (
        lambdaI,
        new volScalarField
        (
            IOobject
            (
                "aLambda_" + Foam::name(lambdaI) ,
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            a_
        )
    );
}

Here, the aLambda_ is initialized by a_, which is initialized in fvDOM constructor as zero value field. The following code will calculate the maximum solid angle omegaMax_:

1
2
3
4
5
6
7
8
// Calculate the maximum solid angle
forAll(IRay_, rayId)
{
    if (omegaMax_ <  IRay_[rayId].omega())
    {
        omegaMax_ = IRay_[rayId].omega();
    }
}

The following code will print Information about what we have done in the discretization:

1
2
3
4
5
6
7
8
9
10
11
Info<< typeName << ": Created " << IRay_.size() << " rays with average "
    << "directions (dAve) and solid angles (omega)" << endl;
Info<< incrIndent;
forAll(IRay_, rayId)
{
    Info<< indent
        << "Ray " << IRay_[rayId].I().name() << ": "
        << "dAve = " << IRay_[rayId].dAve() << ", "
        << "omega = " << IRay_[rayId].omega() << endl;
}
Info<< decrIndent << endl;

Then, we will see how fireFoam call the radiation model. In fireFoam.C, we have #include "YEEqn.H". In YEEqn.H , radiation->correct(); will call the correct() function in the radiation model we choose, which is fvDOM here. But in fvDOM class, there is no function called correct(). But we know, the fvDOM class is a subclass of radiationModel. In class radiationModel, a correct() function exists. After calling the correct() function, a source term radiation->Sh(thermo, he) will be added to the energy equation. Here we will see how we call the correct() function first. in radiationModel.C:

1
2
3
4
5
6
7
8
9
10
11
12
13
void Foam::radiationModel::correct()
{
    if (firstIter_ || (time_.timeIndex() % solverFreq_ == 0))
    {
        calculate();
        firstIter_ = false;
    }

    if (!soot_.empty())
    {
        soot_->correct();
    }
}

Firstly, calculate() will be called and then soot_->correct(); will be called. As we know, the calculate() function is a pure virtual function. As the polymorphism in C++, it will called the calculate() function of its subclass. Here the subclass is the class fvDOM. So the calculate() function in fvDOM will be called:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
void Foam::radiationModels::fvDOM::calculate()
{
    absorptionEmission_->correct(a_, aLambda_);//before this, a_ and aLambda_ are all zeros

    updateBlackBodyEmission();

    // Set rays converged false
    List<bool> rayIdConv(nRay_, false);

    scalar maxResidual = 0;
    label radIter = 0;
    do
    {
        Info<< "Radiation solver iter: " << radIter << endl;

        radIter++;
        maxResidual = 0;
        forAll(IRay_, rayI)
        {
            if (!rayIdConv[rayI])
            {
                scalar maxBandResidual = IRay_[rayI].correct();
                maxResidual = max(maxBandResidual, maxResidual);

                if (maxBandResidual < tolerance_)
                {
                    rayIdConv[rayI] = true;
                }
            }
        }

    } while (maxResidual > tolerance_ && radIter < maxIter_);

    updateG();
}

Firstly, a correct function is called. Here, the absorptionEmission_ actually points to greyMeanCombustion which is a subclass of greyMean. And greyMean is a subclass of absorptionEmissionModel. But in both grepMean and greyMeanCombustion have no function called correct. Therefore, the correct function in absorptionEmissionModel will be called:

1
2
3
4
5
6
7
8
9
void Foam::radiationModels::absorptionEmissionModel::correct
(
    volScalarField& a,
    PtrList<volScalarField>& aj
) const
{
    a = this->a(); // !!! this->a() is in class absorptionEmissionModel!!!
    aj[0] =  a;
}

The thing is weird here, only the first element of aLambda_ which corresponds to parameter aj has been initialized by a_ corresponding to parameter a. If only one wave length is considered in the simulation, this will not influence the result. Then we call updateBlackBodyEmission() which is defined in fvDOM.C:

1
2
3
4
5
6
7
void Foam::radiationModels::fvDOM::updateBlackBodyEmission()
{
    for (label j=0; j < nLambda_; j++)
    {
        blackBody_.correct(j, absorptionEmission_->bands(j));
    }
}

The nLambda_ here is one in this simulation. The blackBody_ is a member of the fvDOM class which has a type of blackBodyEmission. It has been initialized in the constructor of fvDOM as blackBody_(nLambda_, T),. This will call the constructor of the blackBodyEmission:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
Foam::radiationModels::blackBodyEmission::blackBodyEmission
(
    const label nLambda,
    const volScalarField& T
)
:
    table_
    (
        emissivePowerTable,
        interpolationTable<scalar>::CLAMP,
        "blackBodyEmissivePower"
    ),
    C1_("C1", dimensionSet(1, 4, 3, 0, 0, 0, 0), 3.7419e-16),
    C2_("C2", dimensionSet(0, 1, 0, 1, 0, 0, 0), 14.388e-6),
    bLambda_(nLambda),
    T_(T)
{
    forAll(bLambda_, lambdaI)
    {
        bLambda_.set
        (
            lambdaI,
            new volScalarField
            (
                IOobject
                (
                    "bLambda_" + Foam::name(lambdaI) ,
                    T.mesh().time().timeName(),
                    T.mesh(),
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                physicoChemical::sigma*pow4(T)

            )
        );
    }
}

Here the bLambda_ which is the black body emission energy field for each wavelength will be initialized as physicoChemical::sigma*pow4(T), . Here is the Stefan–Boltzmann constant. Then, we go back to the updateBlackBodyEmission() function. The correct function of blackBody_ will be called:

1
2
3
4
5
6
7
8
void Foam::radiationModels::blackBodyEmission::correct
(
    const label lambdaI,
    const Vector2D<scalar>& band
)
{
    bLambda_[lambdaI] = EbDeltaLambdaT(T_, band);
}

In this correct function, the member bLambda_ of blackBody_ will be assigned with EbDeltaLambdaT(T_, band) at wave length lambdaI. The EbDeltaLambdaT function is defined as:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Foam::tmp<Foam::volScalarField>
Foam::radiationModels::blackBodyEmission::EbDeltaLambdaT
(
    const volScalarField& T,
    const Vector2D<scalar>& band
) const
{
    tmp<volScalarField> Eb
    (
        volScalarField::New
        (
            "Eb",
            physicoChemical::sigma*pow4(T)
        )
    );

    if (band != Vector2D<scalar>::one)
    {
     //......
    }

    return Eb;
}

Because here the band we transfer trace from absorptionEmission_->bands(j) in updateBlackBodyEmission function, it means we have to see what band(j) is. In both greyMean and greyMeanCombustion class (greyMean is a subclass of absorptionEmissionModel and greyMeanCombustion is a subclass of greyMean), there is no function called bands, but in absorptionEmissionModel, this bands function exists. Therefore, the bands function called here is the one in the absorptionEmissionModel:

1
2
3
4
5
const Foam::Vector2D<Foam::scalar>&
Foam::radiationModels::absorptionEmissionModel::bands(const label n) const
{
    return Vector2D<scalar>::one;
}

As the returned value is Vector2D<scalar>::one, the if{} will no be called in the EbDeltaLambdaT. Therefore the return value from the EbDeltaLambdaT function here is the Eb, which is equal to physicoChemical::sigma*pow4(T) (). One more thing to mentioned here is what is returned from Vector2D<scalar>::one:

1
Vector2D<scalar>::one // returns (1 1)

Let’s go back to the calculate function of fvDOM. A list of bool with a size of number of rays is defined for convergence check of each radiation ray. The maximum residual maxResidual and a counter for iterations radIter are also defined. These 2 values will be compared with the tolerance and maxIter we set in the radationPrperties file to see if the calculation is converged.

1
2
3
4
// Set rays converged false
List<bool> rayIdConv(nRay_, false);
scalar maxResidual = 0;
label radIter = 0;

In the following do while loop, a condition maxResidual > tolerance_ && radIter < maxIter_ has been set. If the calculation is converged or it reaches the user defined maximum iterator times. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
do
{
    Info<< "Radiation solver iter: " << radIter << endl;

    radIter++;
    maxResidual = 0;
    forAll(IRay_, rayI)
    {
        if (!rayIdConv[rayI])
        {
            scalar maxBandResidual = IRay_[rayI].correct();
            maxResidual = max(maxBandResidual, maxResidual);

            if (maxBandResidual < tolerance_)
            {
                rayIdConv[rayI] = true;
            }
        }
    }

} while (maxResidual > tolerance_ && radIter < maxIter_);

Inside the loop, for each radiation ray IRay_[rayI], the correct function will be called. IRay_ is a pointer list with a type of radiativeIntensityRay. Here, IRay_[rayI] dereference the rayI‘s pointer and call the correct() function of the radiativeIntensityRay class:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Foam::scalar Foam::radiationModels::radiativeIntensityRay::correct()
{
    // Reset boundary heat flux to zero
    qr_.boundaryFieldRef() = 0.0;

    scalar maxResidual = -great;

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

    forAll(ILambda_, lambdaI)
    {
        const volScalarField& k = dom_.aLambda(lambdaI);

        fvScalarMatrix IiEq
        (
            fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
          + fvm::Sp(k*omega_, ILambda_[lambdaI])
        ==
            1.0/constant::mathematical::pi*omega_
           *(
                // Remove aDisp from k
                (k - absorptionEmission_.aDisp(lambdaI))
               *blackBody_.bLambda(lambdaI)

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

        IiEq.relax();

        const solverPerformance ILambdaSol = solve(IiEq, "Ii");

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

        maxResidual = max(initialRes, maxResidual);
    }

    return maxResidual;
}

The mesh_.Sf() returns the surface area vector of the cell faces. Here, const surfaceScalarField Ji(dAve_ & mesh_.Sf()); calculate the vector dot product of dAve_ and mesh_.Sf() as a result of Ji. Then, solve the RTE for radiation ray in direction rayI and wave length lambdaI. Here, we will translate the code of the equation to mathematic equation:

, where

It can be found that the particle scattering factor is not considered in the equation. It will be added to the code later. In the correct() function, the maximum residual of the calculation will also be returned in the end to see if the residual convergence is reached. It should be noted that in the calculate function of the fvDOM class, such a loop will select only the unconverged ray equation to solve. If the ray we are solving is converged, in the rayIdConv, its value will become true and will not be solved again. After this do while loop, the last function to call is the updateG():

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void Foam::radiationModels::fvDOM::updateG()
{
    G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0);
    qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0);
    qem_ = dimensionedScalar(dimMass/pow3(dimTime), 0);
    qin_ = dimensionedScalar(dimMass/pow3(dimTime), 0);

    forAll(IRay_, rayI)
    {
        IRay_[rayI].addIntensity();
        G_ += IRay_[rayI].I()*IRay_[rayI].omega();
        qr_.boundaryFieldRef() += IRay_[rayI].qr().boundaryField();
        qem_.boundaryFieldRef() += IRay_[rayI].qem().boundaryField();
        qin_.boundaryFieldRef() += IRay_[rayI].qin().boundaryField();
    }
}

The G_ is the Incident radiation, which is as same as the G_ in P1 model. Before assigning value to the variables, all of them are initialized to zero with certain dimension. Firstly, addIntensity() is called:

1
2
3
4
5
6
7
8
9
void Foam::radiationModels::radiativeIntensityRay::addIntensity()
{
    I_ = dimensionedScalar(dimMass/pow3(dimTime), 0);

    forAll(ILambda_, lambdaI)
    {
        I_ += ILambda_[lambdaI];
    }
}

Here, the total radiative intensity I_ will be calculated as a summation over all wave length. G_ is calculated as a summation of IRay_[rayI].I()*IRay_[rayI].omega() where I() is defined as: 

1
2
3
4
5
inline const Foam::volScalarField&
Foam::radiationModels::radiativeIntensityRay::I() const
{
    return I_;
}

and omega() is defined as:

1
2
3
4
inline Foam::scalar Foam::radiationModels::radiativeIntensityRay::omega() const
{
    return omega_;
}

Therefore G_ is actually calculated as:

The qr_, qin, andqem are total radiative heat flux on boundary, incident radiative heat flux on boundary, and emitted radiative heat flux on boundary. These values will be updated by a summation over different rays on the boundaries. Till here, the calculate() function in Foam::radiationModel::correct() is finished. Then the radiation->Sh(thermo, he) will be added to the energy equation. As same as we introduced in P1 model previously, the Sh function only exits in the radiationModel class:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Foam::tmp<Foam::fvScalarMatrix> Foam::radiationModel::Sh
(
    const basicThermo& thermo,
    const volScalarField& he
) const
{
    const volScalarField Cpv(thermo.Cpv());
    const volScalarField T3(pow3(T_));

    return
    (
        Ru()
      - fvm::Sp(4.0*Rp()*T3/Cpv, he)
      - Rp()*T3*(T_ - 4.0*he/Cpv)
    );
}

The only difference is how we calculate the Ru() and Rp(). In fvDOM, Ru() is calculated as:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
Foam::tmp<Foam::volScalarField> Foam::radiationModels::fvDOM::Rp() const
{
    // Construct using contribution from first frequency band
    tmp<volScalarField> tRp
    (
        volScalarField::New
        (
            "Rp",
            4
           *physicoChemical::sigma
           *(aLambda_[0] - absorptionEmission_->aDisp(0)())
           *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(0))
        )
    );

    volScalarField& Rp=tRp.ref();

    // Add contributions over remaining frequency bands
    for (label j=1; j < nLambda_; j++)
    {
        Rp +=
        (
            4
           *physicoChemical::sigma
           *(aLambda_[j] - absorptionEmission_->aDisp(j)())
           *blackBody_.deltaLambdaT(T_, absorptionEmission_->bands(j))
        );
    }

    return tRp;
}

The function deltaLambdaT is defined as in class blackBodyEmission:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
Foam::tmp<Foam::volScalarField>
Foam::radiationModels::blackBodyEmission::deltaLambdaT
(
    const volScalarField& T,
    const Vector2D<scalar>& band
) const
{
    tmp<volScalarField> deltaLambdaT
    (
        volScalarField::New
        (
            "deltaLambdaT",
            T.mesh(),
            dimensionedScalar(dimless, 1.0)
        )
    );

    if (band != Vector2D<scalar>::one)
    {
        scalarField& deltaLambdaTf = deltaLambdaT.ref();

        forAll(T, i)
        {
            deltaLambdaTf[i] = fLambdaT(band[1]*T[i]) - fLambdaT(band[0]*T[i]);
        }
    }

    return deltaLambdaT;
}

We only have one band, there for deltaLambdaT is one. Rp() is calculated as the following equation:

The definition of Ru is:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
Foam::radiationModels::fvDOM::Ru() const
{
    tmp<DimensionedField<scalar, volMesh>> tRu
    (
        new DimensionedField<scalar, volMesh>
        (
            IOobject
            (
                "Ru",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh_,
            dimensionedScalar(dimensionSet(1, -1, -3, 0, 0), 0)
        )
    );

    DimensionedField<scalar, volMesh>& Ru=tRu.ref();

    // Sum contributions over all frequency bands
    for (label j=0; j < nLambda_; j++)
    {
        // Compute total incident radiation within frequency band
        tmp<DimensionedField<scalar, volMesh>> Gj
        (
            IRay_[0].ILambda(j)()*IRay_[0].omega()
        );

        for (label rayI=1; rayI < nRay_; rayI++)
        {
            Gj.ref() += IRay_[rayI].ILambda(j)()*IRay_[rayI].omega();
        }

        Ru += (aLambda_[j]() - absorptionEmission_->aDisp(j)()())*Gj
             - absorptionEmission_->ECont(j)()();
    }

    return tRu;
}

The code follows equation:

Then we can substitute the Ru() and Rp() into the Sh and add them to the energy equation as a source term.