P1 Model

To understand the how P1 model is used and implemented in coalChemistryFoam, we use tutorial simplifiedSiwek as an example.

P1 model’s euqtion:

In the code, it looks like this:

1
2
3
4
5
6
7
solve
(
    fvm::laplacian(gamma, G_)
  - fvm::Sp(a_, G_)
 ==
  - 4.0*(e_*physicoChemical::sigma*pow4(T_) ) - E_
);

which is defined in Foam::radiationModels::P1::calculate in P1.C. Here, equals gamma (the diffusivity), equals a_ (the absorption coefficient ), equals G_, equals e_ (the emission coefficient), equals pow4(T_), equals E_ (the emission contribution). is the Stefan-Boltzmann constant. The definition of follows:

where is a infinitely value to make sure denominator is not equal to zero in the simulation.

Before starting, the radiationProperties file used in tutorial simplifiedSiwek will be listed here:

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
solverFreq      1;

radiationModel  P1;

absorptionEmissionModel binary;

binaryCoeffs
{
    model1
    {
        absorptionEmissionModel constant;
        constantCoeffs
        {
            absorptivity    0.5; // a_
            emissivity      0.5; // e_
            E               0; // user defined radiation source (the emission contribution)
            boyao           zixin;
            boyao1          gouzi;
        }
    }
    model2
    {
        absorptionEmissionModel cloud;
        cloudCoeffs
        {
            cloudNames
            (
                coalCloud1
                limestoneCloud1
            );
        }
    }
}

scatterModel    cloud;

cloudCoeffs
{
    cloudNames
    (
        coalCloud1
        limestoneCloud1
    );
}

In coalChemistryFoam.C

createFields.H is included. In createFields.H, createRadiationModel.H is included. In createRadiationModel.H:

1
autoPtr<radiationModel> radiation(radiationModel::New(thermo.T()));

a object called radiation is constructed, which has type of a pointer to class type radiationModel. autoPtr is similar to a smart pointer in C++. This radiationModel object is called in EEqn.H radiation->correct(); Of course, in calss radiationModel, there is a member function called virtual void correct() which is a virtual function. In its correct() function, it calls calculate() which in this class is a pure virtual function. It means this function will be called in its sub classes. Of course, P1::calculate() will be called in this case. In EEqn.H, radiation->Sh(thermo, he) is also called.

In radiationModel folder, there is a file called radiationModelNew.C. In this file, Foam::autoPtr<Foam::radiationModel> Foam::radiationModel::New(const volScalarField& T) is defined.

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
Foam::autoPtr<Foam::radiationModel> Foam::radiationModel::New
(
    const volScalarField& T
)
{
    IOobject radIO
    (
        "radiationProperties", //this indicates which file in constant folder it is in
        T.time().constant(),
        T.mesh(),
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE,
        false
    );

    word modelType("none"); 
    if (radIO.typeHeaderOk<IOdictionary>(false))
    {
        IOdictionary(radIO).lookup("radiationModel") >> modelType; //here you search for the model e.g. P1 model in the run folder
    }
    else
    {
        Info<< "Radiation model not active: radiationProperties not found"
            << endl;
    }

    Info<< "Selecting radiationModel " << modelType << endl; // print the model you have chosed

    TConstructorTable::iterator cstrIter =
        TConstructorTablePtr_->find(modelType);
        
    //***code below can be inserted here to see which models are in the table***
    
    // for(auto it=dictionaryConstructorTablePtr_->begin();
    //     it!= dictionaryConstructorTablePtr_->end();it++ )
    //     {
    //       cout << it.key() << '\n';
    //     }

    if (cstrIter == TConstructorTablePtr_->end())
    {
        FatalErrorInFunction
            << "Unknown radiationModel type "
            << modelType << nl << nl
            << "Valid radiationModel types are:" << nl
            << TConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }

    return autoPtr<radiationModel>(cstrIter()(T));
}

The modelType here has a type of word and is initialised as "P1". This is done by IOdictionary(radIO).lookup("radiationModel") >> modelType. Finally, the constructor of the radiationModel will be called in radiationModel.C. This is done by the last line of the Foam::autoPtr<Foam::radiationModel> Foam::radiationModel::New. Here a smart pointer with a type radiationModel is returned. This pointer is initialized by cstrIter()(T). Here, cstrIter has a type of TConstructorTable::iterator where TConstructorTable is exactly the HashTable type, defined by run time selection Macro. Therefore, cstrIter has a type of HashTable<dictionaryConstructorPtr, word, string::hash>::iterator. The iterator class is defined inside the HashTable class, it has an operator() defined in the class body which is called here (cstrIter()(T)). Here, the cstrIter() actually represent the P1 model. Therefore, return autoPtr<radiationModel>(cstrIter()(T)); will call the constructor of radiationModels::P1 class, which is a subclass of Foam::radiationModel:

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
Foam::radiationModels::P1::P1(const volScalarField& T)
:
    radiationModel(typeName, T),
    G_
    (
        IOobject
        (
            "G",
            mesh_.time().timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),
    qr_
    (
        IOobject
        (
            "qr",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar(dimMass/pow3(dimTime), 0)
    ),
    a_
    (
        IOobject
        (
            "a",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar(dimless/dimLength, 0)
    ),
    e_
    (
        IOobject
        (
            "e",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar(dimless/dimLength, 0)
    ),
    E_
    (
        IOobject
        (
            "E",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar(dimMass/dimLength/pow3(dimTime), 0)
    )
{}

In the constructor of radiationModels::P1 class, it will first initialize its parent class radiationModel by calling its constructor(shown below). Then, its members (G_, qr_, a_, e_, E_) will be initialized. The definitions of those variables are shown below:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
//- Incident radiation / [W/m^2]
volScalarField G_;
 
//- Total radiative heat flux [W/m^2]
volScalarField qr_;
 
//- Absorption coefficient
volScalarField a_;
 
//- Emission coefficient
volScalarField e_;
 
//- Emission contribution
volScalarField E_;

The constructor of the Foam::radiationModel is shown 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 use of IOobject class has been explained here. There is a initialise() function called in the constructor which is also defined in radiationModel.C.

1
2
3
4
5
6
7
8
9
10
11
12
13
void Foam::radiationModel::initialise()
{
    solverFreq_ = max(1, lookupOrDefault<label>("solverFreq", 1));

    absorptionEmission_.reset
    (
        radiationModels::absorptionEmissionModel::New(*this, mesh_).ptr()
    );

    scatter_.reset(radiationModels::scatterModel::New(*this, mesh_).ptr());

    soot_.reset(radiationModels::sootModel::New(*this, mesh_).ptr());
}

Then we go back to initialise(). It shoud be noted here, the class radiationModel is a subclass of IOdictionary which is a subclass of dictionary. Therefore when we called e.g. radiationModels::absorptionEmissionModel::New(*this, mesh_) which is defined in absorptionEmissionModelNew.C:

1
2
3
4
5
6
 Foam::autoPtr<Foam::radiationModels::absorptionEmissionModel>
Foam::radiationModels::absorptionEmissionModel::New
(
    const dictionary& dict,
    const fvMesh& mesh
)

the *this which is a radiationModel object is substituted. absorptionEmission_, scatter_, scatter_ are defined in radiationModel.H:

1
2
3
4
5
6
7
8
autoPtr<radiationModels::absorptionEmissionModel>
    absorptionEmission_;

//- Scatter model
autoPtr<radiationModels::scatterModel> scatter_;

//- Soot model
autoPtr<radiationModels::sootModel> soot_;

In the initialise() function, the three autoPtr smart pointers will call each respected constructors by calling the reset function. The first one:

1
2
3
4
absorptionEmission_.reset
(
    radiationModels::absorptionEmissionModel::New(*this, mesh_).ptr()
);

will call the absorptionEmissionModel::New, which is quite similar to Foam::radiationModel::New:

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
Foam::autoPtr<Foam::radiationModels::absorptionEmissionModel>
Foam::radiationModels::absorptionEmissionModel::New
(
    const dictionary& dict,
    const fvMesh& mesh
)
{
    const word modelType(dict.lookup("absorptionEmissionModel"));

    Info<< "Selecting absorptionEmissionModel " << modelType << endl;

    dictionaryConstructorTable::iterator cstrIter =
        dictionaryConstructorTablePtr_->find(modelType);

    if (cstrIter == dictionaryConstructorTablePtr_->end())
    {
        FatalErrorInFunction
            << "Unknown absorptionEmissionModel type "
            << modelType << nl << nl
            << "Valid absorptionEmissionModel types are :" << nl
            << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError);
    }

    return autoPtr<absorptionEmissionModel>(cstrIter()(dict, mesh));
}

Here in the radiationProperties file, const word modelType(dict.lookup("absorptionEmissionModel")); will result in binary. Therefore, return autoPtr<absorptionEmissionModel>(cstrIter()(dict, mesh)); will call the constructor of class binary in binary.C as binary is a subclass of absorptionEmissionModel:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Foam::radiationModels::absorptionEmissionModels::binary::binary
(
    const dictionary& dict,
    const fvMesh& mesh
)
:
    absorptionEmissionModel(dict, mesh),
    coeffsDict_(dict.optionalSubDict(typeName + "Coeffs")),
    model1_
    (
        absorptionEmissionModel::New(coeffsDict_.subDict("model1"), mesh)
    ),
    model2_
    (
        absorptionEmissionModel::New(coeffsDict_.subDict("model2"), mesh)
    )
{}

Here the absorptionEmissionModel::New will be called again. The use of subDict and optionalSubDict has been explained here. The typeName here is binary, which means typeName + "Coeffs" is binaryCoeffs. Therefore, coeffsDict_ will be:

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
{
    model1
    {
        absorptionEmissionModel constant;
        constantCoeffs
        {
            absorptivity    0.5;
            emissivity      0.5;
            E               0;
            boyao           zixin;
            boyao1          gouzi;
        }
    }
    model2
    {
        absorptionEmissionModel cloud;
        cloudCoeffs
        {
            cloudNames
            (
                coalCloud1
                limestoneCloud1
            );
        }
    }
}

So, coeffsDict_.subDict("model1") in this case will be a dictionary of:

1
2
3
4
5
6
7
8
9
10
11
{
    absorptionEmissionModel constant;
    constantCoeffs
    {
        absorptivity    0.5;
        emissivity      0.5;
        E               0;
        boyao           zixin;
        boyao1          gouzi;
    }
}

Vice versa the coeffsDict_.subDict("model2"). Let’s go forward. The binary constructor will also initialize its members model1_ and model2_, which are defined in binary.H:

1
2
3
4
5
//- First absorption model
autoPtr<absorptionEmissionModel> model1_;
 
//- Second absorption model
autoPtr<absorptionEmissionModel> model2_;

These 2 members have a type of pointer to absorptionEmissionModel. When we initialize it, it will call the absorptionEmissionModel::New function again. Now we will go through this function again to make everything clear. When we call const word modelType(dict.lookup("absorptionEmissionModel")); in absorptionEmissionModel::New function, it will return a constant to modelType. Then, return autoPtr<absorptionEmissionModel>(cstrIter()(dict, mesh)); will call the constructor of class constant in constantAbsorptionEmission.C:

1
2
3
4
5
6
7
8
9
10
11
12
13
Foam::radiationModels::absorptionEmissionModels::constant::constant
(
    const dictionary& dict,
    const fvMesh& mesh
)
:
    absorptionEmissionModel(dict, mesh),
    coeffsDict_(dict.optionalSubDict(typeName + "Coeffs")),
    a_("absorptivity", dimless/dimLength, coeffsDict_),
    e_("emissivity", dimless/dimLength, coeffsDict_),
    E_("E", dimMass/dimLength/pow3(dimTime), coeffsDict_)
{}

Here, we will first call the constructor of absorptionEmissionModel. This is the first time we call the constructor of the base class absorptionEmissionModel:

1
2
3
4
5
6
7
8
9
Foam::radiationModels::absorptionEmissionModel::absorptionEmissionModel
(
    const dictionary& dict,
    const fvMesh& mesh
)
:
    dict_(dict),
    mesh_(mesh)
{}

This constructor will initialize two protected members declared in absorptionEmissionModel.H:

1
2
3
4
5
6
7
8
9
protected:

    // Protected data

        //- Radiation model dictionary
        const dictionary dict_;

        //- Reference to the fvMesh
        const fvMesh& mesh_;

Those protected members can be used in its subclass e.g. constant. Then the coeffsDict_ will be initialized as coeffsDict_(dict.optionalSubDict(typeName + "Coeffs")), where typeName here is constant. Therefore coeffsDict_ is:

1
2
3
4
5
6
7
8
constantCoeffs
{
    absorptivity    0.5;
    emissivity      0.5;
    E               0;
    boyao           zixin;
    boyao1          gouzi;
}

Then, a_, e_ E_ will be initialized with the values in absorptivity, emissivity, E. For the initialization of model2_, there is only a little bit difference from that of model1_. Because the model2_ is initialized with type cloud. It is defined in a totally different folder (src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission). This class which is called cloud is also a subclass of absorptionEmissionModel defined in cloudAbsorptionEmission.H. In cloudAbsorptionEmission.C:

1
2
3
4
5
6
7
8
9
10
11
Foam::radiationModels::absorptionEmissionModels::cloud::cloud
(
    const dictionary& dict,
    const fvMesh& mesh
)
:
    absorptionEmissionModel(dict, mesh),
    coeffsDict_(dict.subDict(typeName + "Coeffs")),
    cloudNames_(coeffsDict_.lookup("cloudNames"))
{}

Here, no further explanation is needed for how the constructor is working. Due to its similarity to the constructor of class constant. For scatter model, it will be similar to absorptionEmissionModel, we won’t go through it again.

In coalChemistryFoam.C we include EEqn.H. In EEqn.H, the radiation object which was constructed will be used. It is used as radiation->Sh(thermo, he) and radiation->correct(). The radiation->correct() calls void Foam::radiationModels::P1::calculate(). One thing here we have to be pretty clear about is where the pointer absorptionEmission_ is pointing to !!!!!!! As we do the .reset in the initialise() function, we call the radiationModels::absorptionEmissionModel::New(*this, mesh_) which return a pointer points to binary class. Therefore, here the absorptionEmission_ pointing to the binary class. So I will note here once more, the absorptionEmission_ points to the binary class. Here, I will show the Foam::radiationModels::P1::calculate() of P1.C:

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
void Foam::radiationModels::P1::calculate()
{
    a_ = absorptionEmission_->a();
    e_ = absorptionEmission_->e();
    E_ = absorptionEmission_->E();
    const volScalarField sigmaEff(scatter_->sigmaEff());

    const dimensionedScalar a0 ("a0", a_.dimensions(), rootVSmall); 
    // a infinitely small value to keep the denominator larger than 0

    // Construct diffusion
    const volScalarField gamma
    (
        IOobject
        (
            "gammaRad",
            G_.mesh().time().timeName(),
            G_.mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        1.0/(3.0*a_ + sigmaEff + a0)
    );

    // Solve G transport equation
    solve
    (
        fvm::laplacian(gamma, G_)
      - fvm::Sp(a_, G_)
     ==
      - 4.0*(e_*physicoChemical::sigma*pow4(T_) ) - E_
    );

    volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();

    // Calculate radiative heat flux on boundaries.
    forAll(mesh_.boundaryMesh(), patchi)
    {
        if (!G_.boundaryField()[patchi].coupled())
        {
            qrBf[patchi] =
                -gamma.boundaryField()[patchi]
                *G_.boundaryField()[patchi].snGrad();
        }
    }
}

First, the calculate function initialize the volScalarField members we need e.g. a_, e_, E_, a0, gamma. Then we will solve the P1 equation in solve. I will explain here how a_ is initialized in this calculate() function to represent how these volScalarField are initialized. Here, I will note it once more, the absorptionEmission_ is pointing to binary class, which is a subclass of absorptionEmission class. absorptionEmission_->a() will call function a(). But a() only exist in absorptionEmission class. Therefore we will call this a() in absorptionEmission class first. In absorptionEmissionModel.H:

1
2
//- Absorption coefficient (net)
virtual tmp<volScalarField> a(const label bandI = 0) const;

In absorptionEmissionModel.C:

1
2
3
4
5
Foam::tmp<Foam::volScalarField>
Foam::radiationModels::absorptionEmissionModel::a(const label bandI) const
{
    return aDisp(bandI) + aCont(bandI);
}

In aDisp(bandI) + aCont(bandI), function aDisp and aCont in class binary will be called. function aDisp and aCont are virtual functions firstly defined in absorptionEmissionModel class. As we mentioned, binary is a subclass of absorptionEmissionModel and it also has the aDisp and aCont functions. As absorptionEmission_ is pointing to class binary, the aDisp and aCont function in binary class will be called. I will show aCont function as an example here to see how the function is called. In binary.H:

1
2
//- Absorption coefficient for continuous phase
virtual tmp<volScalarField> aCont(const label bandI = 0) const;

In binary.C:

1
2
3
4
5
6
7
8
Foam::tmp<Foam::volScalarField>
Foam::radiationModels::absorptionEmissionModels::binary::aCont
(
    const label bandI
) const
{
    return model1_->aCont(bandI) + model2_->aCont(bandI);
}

Members model1_ and model2_ are defined here. Both are smart pointers. model1_ points to class constant while model2_ points to class radiationModels::absorptionEmissionModels::cloud. In constantAbsorptionEmission.C:

1
2
3
4
5
6
7
8
9
10
11
12
13
Foam::tmp<Foam::volScalarField>
Foam::radiationModels::absorptionEmissionModels::constant::aCont
(
    const label bandI
) const
{
    return volScalarField::New
    (
        "a",
        mesh_,
        a_   // from input
    );
}

It means aCont will be initialized with the a_ we input in the radiationProperties file. a_ is initialized here (in the constant constructor). But in class Foam::radiationModels::absorptionEmissionModels::cloud , there is no function called aCont. It is reasonable, because in particle field, there shouldn’t be a aCont function. As virtual function’s polymorphism in C++, model2_->aCont will call the aCont() function in its base class absorptionEmissionModel:

1
2
3
4
5
6
7
8
9
10
Foam::tmp<Foam::volScalarField>
Foam::radiationModels::absorptionEmissionModel::aCont(const label bandI) const
{
    return volScalarField::New
    (
        "aCont",
        mesh_,
        dimensionedScalar(dimless/dimLength, 0)
    );
}

As you see, it will be initialized with zeros in this volScalarField.

After all these initialization, we will solve the P1 model’s equation in solve() and return a source term to the energy equation in EEqn.H:

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
{
    volScalarField& he = thermo.he();

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        rho*(U&g)
      + combustion->Qdot()
      + coalParcels.Sh(he)
      + limestoneParcels.Sh(he)
      + radiation->Sh(thermo, he)
      + fvOptions(rho, he)
    );

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    fvOptions.correct(he);

    thermo.correct();
    radiation->correct();

    Info<< "T gas min/max   = " << min(T).value() << ", "
        << max(T).value() << endl;
}

The Sh() is defined in radiationModel, which is not existing in P1 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 calculation of Ru and Rp are defined in P1 as they are virtual functions both in radiationModel and P1 class:

1
2
3
4
5
6
7
8
9
10
11
12
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
Foam::radiationModels::P1::Ru() const
{
    const volScalarField::Internal& G =
        G_();
    const volScalarField::Internal E =
        absorptionEmission_->ECont()()();
    const volScalarField::Internal a =
        absorptionEmission_->aCont()()();

    return a*G - E;
}

and

1
2
3
4
5
6
7
8
Foam::tmp<Foam::volScalarField> Foam::radiationModels::P1::Rp() const
{
    return volScalarField::New
    (
        "Rp",
        4.0*absorptionEmission_->eCont()*physicoChemical::sigma
    );

Here we will explain the model’s equations e.g. how the radiation term acts on the energy transport equation. The energy source term returned to the heat transport equation here is the Sh(). It is calculated as(in the code):

where,

results in,

In ANSYS FLUENT, the form of the equation is actually the same:

It should be mentioned here, is zero. So we don’t need to multiply the total .

radiationModel
radiationModel.C
absorptionEmissionModel
absorptionEmissionModel.C
absorptionEmissionModel.H
P1.C