Devolatilization Model in OpenFOAM

Devolatilization Model in OpenFOAM

We are going to talk about how the devolatilization model used in OpenFOAM. It is defined in the constant/coalCloudProperties. E.g. in the tutorial simplifiedSiwek, which has used the coalChemistryFoam.

First, we find that in coalChemistryFoam.C, it includes the coalCloud.H. In coalCloud.H, it defines:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
typedef ReactingMultiphaseCloud
<
    ReactingCloud
    <
        ThermoCloud
        <
            KinematicCloud
            <
                Cloud
                <
                    coalParcel
                >
            >
        >
    >
> coalCloud;

which is a highly inherited template (see template inheritance). The inner most class inside the angle brackets is called coalParcel which is defined in coalParcel.H:

1
2
3
4
5
6
7
8
9
10
11
12
13
typedef ReactingMultiphaseParcel
<
    ReactingParcel
    <
        ThermoParcel
        <
            KinematicParcel
            <
                particle
            >
        >
    >
> coalParcel;

, which is also a template inheritance class. The inner most part is called particle which is a class defined in particle.H.

Let’s go back to the coalChemistyFoam. In coalChemistyFoam, it includes a header called createFields.H, which includes createClouds.H:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Info<< "\nConstructing coal cloud" << endl;
coalCloud coalParcels
(
    "coalCloud1",
    rho,
    U,
    g,
    slgThermo
);

Info<< "\nConstructing limestone cloud" << endl;
basicThermoCloud limestoneParcels
(
    "limestoneCloud1",
    
    rho,
    U,
    g,
    slgThermo
);

In the header above, the coalCloud class is used by defining an object called coalParcels. This object will call the constructor of ReactingMultiphaseCloud (Remember here, it can only call the constructor of outer most class of the inherited template) in ReactingMultiphaseCloud.H:

1
2
3
4
5
6
7
8
9
ReactingMultiphaseCloud
(
    const word& cloudName,
    const volScalarField& rho,
    const volVectorField& U,
    const dimensionedVector& g,
    const SLGThermo& thermo,
    bool readFields = true
);

and in ReactingMultiphaseCloud.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
template<class CloudType>
Foam::ReactingMultiphaseCloud<CloudType>::ReactingMultiphaseCloud
(
    const word& cloudName,
    const volScalarField& rho,
    const volVectorField& U,
    const dimensionedVector& g,
    const SLGThermo& thermo,
    bool readFields
)
:
    CloudType(cloudName, rho, U, g, thermo, false),
    reactingMultiphaseCloud(),
    cloudCopyPtr_(nullptr),
    constProps_(this->particleProperties()),
    devolatilisationModel_(nullptr),
    surfaceReactionModel_(nullptr),
    dMassDevolatilisation_(0.0),
    dMassSurfaceReaction_(0.0)
{
    if (this->solution().active())
    {
        setModels();

        if (readFields)
        {
            parcelType::readFields(*this, this->composition());
            this->deleteLostParticles();
        }
    }

    if (this->solution().resetSourcesOnStartup())
    {
        resetSourceTerms();
    }
}

In the constructor of the ReactingMultiphaseCloud, it will initialize the next constructor in one layer inner of the angle bracket. This is why OpenFOAM use this template inheritance method. The outer class can not only inherit the inner class functions but also its type. If you check each layer of the template, you will find out its relationship to each class layer.

The Devolatilisation model object devolatilisationModel_ is declared in ReactingMultiphaseCloud.H as a smart pointer (autoPtr):

1
2
3
4
5
6
7
8
// References to the cloud sub-models
 
    //- Devolatilisation model
    autoPtr
    <
        DevolatilisationModel<ReactingMultiphaseCloud<CloudType>>
    >
    devolatilisationModel_;

One question remained here, which I didn’t figure out is why the Devolatilization class is forward declared:

1
2
template<class CloudType>
class DevolatilisationModel;

I would prefer to include the header file here instead of the forward declaration. In the constructor of the ReactingMultiphaseCloud, the devolatilisationModel_ is initialized as a nullptr. In the constructor body, the function setModels() will be called. 

1
2
3
4
5
6
7
8
9
10
11
12
13
template<class CloudType>
void Foam::ReactingMultiphaseCloud<CloudType>::setModels()
{
    devolatilisationModel_.reset
    (
        DevolatilisationModel<ReactingMultiphaseCloud<CloudType>>::New
        (
            this->subModelProperties(),
            *this
        ).ptr()
    );
    ...
}

In this function, the devolatilizationModel_ , which is a smart pointer will be reset, (the reset is a function used in the shared pointer class, therefore it can be called using .) with a new pointer. The DevolatilisationModel<ReactingMultiphaseCloud<CloudType>>::New is defined in DevolatilisationModelNew.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
template<class CloudType>
Foam::autoPtr<Foam::DevolatilisationModel<CloudType>>
Foam::DevolatilisationModel<CloudType>::New
(
    const dictionary& dict,
    CloudType& owner
)
{
    const word modelType(dict.lookup("devolatilisationModel"));

    Info<< "Selecting devolatilisation model " << modelType << endl;

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

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

    return autoPtr<DevolatilisationModel<CloudType>>(cstrIter()(dict, owner));
}

This kind of New function is very common in OpenFOAM, we have discussed about this in previous posts e.g. P1 Model. In below, we will study how the two parameters are implemented in this New fuction.

The first input parameter in the New function is the this->subModelProperties(),. The subModelProperties() is defined in KinematicCloudI.H:

1
2
3
4
5
6
7
template<class CloudType>
inline const Foam::dictionary&
Foam::KinematicCloud<CloudType>::subModelProperties() const
{
    return subModelProperties_;
}

The reason why the subModelProperties() can be called is because the KinematicCloud is a also inside the template inheritance bracket. Don’t forget that the CloudType here is the ReactingCloud, this can be found in the definition of the inherited template defined above. The subModelProperties_ is defined in KinematicCloud.H:

1
2
//- Sub-models dictionary
const dictionary subModelProperties_;

It should be noted here the IOdictionary is a subclass of the dictionary class. In the constructor of the IOdictionary, the IOobjective can be an input parameter:

1
2
//- Construct given an IOobject
IOdictionary(const IOobject&);

In the constructor of KinematicCloud, the subModelProperties_ is initialized as:

1
2
3
4
subModelProperties_
(
    particleProperties_.subOrEmptyDict("subModels", solution_.active())
),

Then we have to see how particleProperties_ is initialized:

1
2
3
4
5
6
7
8
9
10
11
particleProperties_
(
    IOobject
    (
        cloudName + "Properties",
        rho.mesh().time().constant(),
        rho.mesh(),
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
),

The particleProperties_ is defined in KinematicCloud.H:

1
2
//- Dictionary of particle properties
IOdictionary particleProperties_;

The cloudName here is coalCloud1 as we initialized the coalParcels. Therefore, this particleProperties will find the file coalCloud1Properties in the constant folder. The subOrEmptyDict is defined in dictionary.H:

1
2
3
4
5
6
7
//- Find and return a sub-dictionary as a copy, or
//  return an empty dictionary if the sub-dictionary does not exist
dictionary subOrEmptyDict
(
    const word&,
    const bool mustRead = false
) const;

This subOrEmptyDict function will return the sub-directory, which is the subModels in this case, of the file coalCloud1Properties in the constant folder.

Then we can go back to DevolatilisationModel<CloudType>::New, which has 2 input parameters with types of dictionary and CloudType. It will first lookup the devolatilization model in the subModels:

1
const word modelType(dict.lookup("devolatilisationModel"));

Then find this model in the prepared hash table (use the runtime selection mechanism):

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

If the model is successfully found in the prepared hash table, this new function will return a pointer and the corresponding constructor of the devolatilization model will be called:

1
return autoPtr<DevolatilisationModel<CloudType>>(cstrIter()(dict, owner));

Here in the tutorial, the devolatilization model used is called constantRateDevolatilisation. Therefore the construcXDMDgxlf001!!tor of constantRateDevolatilisation, which is defined in ConstantRateDevolatilisation.C 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
template<class CloudType>
Foam::ConstantRateDevolatilisation<CloudType>::ConstantRateDevolatilisation
(
    const dictionary& dict,
    CloudType& owner
)
:
    DevolatilisationModel<CloudType>(dict, owner, typeName), //typeName is constantRateDevolatilisation
    volatileData_(this->coeffDict().lookup("volatileData")),
    YVolatile0_(volatileData_.size()),
    volatileToGasMap_(volatileData_.size()),
    residualCoeff_(readScalar(this->coeffDict().lookup("residualCoeff")))
{
    if (volatileData_.empty())
    {
        WarningInFunction
            << "Devolatilisation model selected, but no volatiles defined"
            << nl << endl;
    }
    else
    {
        Info<< "Participating volatile species:" << endl;

        // Determine mapping between active volatiles and cloud gas components
        const label idGas = owner.composition().idGas();
        const scalar YGasTot = owner.composition().YMixture0()[idGas];
        const scalarField& YGas = owner.composition().Y0(idGas);
        forAll(volatileData_, i)
        {
            const word& specieName = volatileData_[i].first();
            const label id = owner.composition().localId(idGas, specieName);
            volatileToGasMap_[i] = id;
            YVolatile0_[i] = YGasTot*YGas[id];

            Info<< "    " << specieName << ": particle mass fraction = "
                << YVolatile0_[i] << endl;
        }
    }
}

Because the ConstantRateDevolatilisation class inherits the DevolatilisationModel class whose default constructor does not exist, this constructor needs to be initialized first. The constructor of the DevolatilisationModel is:

1
2
3
4
5
6
7
8
9
10
11
template<class CloudType>
Foam::DevolatilisationModel<CloudType>::DevolatilisationModel
(
    const dictionary& dict,
    CloudType& owner,
    const word& type // type is constantRateDevolatilisation
)
:
    CloudSubModelBase<CloudType>(owner, dict, typeName, type), // typeName is devolatilisationModel 
    dMass_(0.0)
{}

Because DevolatilisationModel inherits from the CloudSubModelBase, the constructor of the CloudSubModelBase needs to be initialized. Here the typeName is the devolatilisationModel as defined in DevolatilisationModel.H. The constructor of the CloudSubModelBase is:

1
2
3
4
5
6
7
8
9
//- Construct from owner cloud without name
CloudSubModelBase
(
    CloudType& owner,
    const dictionary& dict,
    const word& baseName,
    const word& modelType,
    const word& dictExt = "Coeffs"
);

and

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
template<class CloudType>
Foam::CloudSubModelBase<CloudType>::CloudSubModelBase
(
    CloudType& owner,
    const dictionary& dict,
    const word& baseName,
    const word& modelType,
    const word& dictExt
)
:
    subModelBase
    (
        owner.outputProperties(),
        dict,
        baseName,
        modelType,
        dictExt
    ),
    owner_(owner)
{}

The subModelBase constructor:

1
2
3
4
5
6
7
8
9
//- Construct from components without name
subModelBase
(
    dictionary& properties,
    const dictionary& dict,
    const word& baseName,
    const word& modelType,
    const word& dictExt = "Coeffs"
);

and

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Foam::subModelBase::subModelBase
(
    dictionary& properties,
    const dictionary& dict,
    const word& baseName,
    const word& modelType,
    const word& dictExt
)
:
    modelName_(word::null),
    properties_(properties),
    dict_(dict),
    baseName_(baseName),
    modelType_(modelType),
    coeffDict_(dict.subDict(modelType + dictExt))
{}

Till here, all the constructor initialization is finished. All those initialization is due to the inheritance of each class to its parent class and in their parent class, there is not default constructor. If the default constructor of the parent class is not existed, it is necessary to initialize the parent class in the constructor. Then we will go to these four initializations:

1
2
3
4
volatileData_(this->coeffDict().lookup("volatileData")),
YVolatile0_(volatileData_.size()),
volatileToGasMap_(volatileData_.size()),
residualCoeff_(readScalar(this->coeffDict().lookup("residualCoeff")))

They are defined in ConstantRateDevolatilisation.H:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// Model constants
 
    //- List of volatile data - (name A0)
    List<Tuple2<word, scalar>> volatileData_;
 
    //- List of initial volatile mass fractions
    List<scalar> YVolatile0_;
 
    //- Mapping between local and cloud gaseous species
    List<label> volatileToGasMap_;
 
    //- Volatile residual coefficient (0-1)
    //  When the fraction of volatiles are depleted below this
    //  threshold, combustion can occur
    const scalar residualCoeff_;
can store two objects of different types, which is defined in ```Tuple2.H```. In the initialization of ```volatileData_```, the ```coeffDict()```, which is defined in ```subModelBase.C``` is called:
1
2
3
4
5
6

```cpp
 const Foam::dictionary& Foam::subModelBase::coeffDict() const
 {
     return coeffDict_;
 }

The coeffDict_ defined in subModelBase.H has a type of dictionary. It is initialized in the constructor of subModelBase. The dict here is the subModels dictionary in the constant folder coalCloud1Properties. The modelType (type) here is the constantRateDevolatilisation. The modelType + dictExt is constantRateDevolatilisationCoeffs. In the simplifiedSiwek tutorial, the dict.subDict(modelType + dictExt) will be:

1
2
3
4
5
6
7
8
9
{
    volatileData
    (
        (CH4            12)
        (H2             12)
        (CO2            12)
    );
    residualCoeff   0.001;
}

The this->coeffDict().lookup("volatileData") will be:

1
2
3
(CH4            12)
(H2             12)
(CO2            12)

, which are three tuples. As volatileData_ is a list, its size here will be 3. YVolatile0_ and volatileToGasMap_ are initialized by size. This can be found here. The initialization of residualCoeff_ is similar to that of volatileData_, it will not be discussed here again. Then we dig into how they are used:

1
2
3
const label idGas = owner.composition().idGas();
const scalar YGasTot = owner.composition().YMixture0()[idGas];
const scalarField& YGas = owner.composition().Y0(idGas);

The owner here is the *this of the ReactingMultiphaseCloud class. The composition() function is defined in ReactingCloud.H:

1
2
3
//- Return const access to reacting composition model
inline const CompositionModel<ReactingCloud<CloudType>>&
    composition() const;

and in ReactingCloudI.H:

1
2
3
4
5
6
template<class CloudType>
inline const Foam::CompositionModel<Foam::ReactingCloud<CloudType>>&
Foam::ReactingCloud<CloudType>::composition() const
{
    return compositionModel_;
}

The compositionModel is defined in ReactingCloud.H:

1
2
3
//- Reacting composition model
autoPtr<CompositionModel<ReactingCloud<CloudType>>>
    compositionModel_;

which is initialized in the constructor of ReactingCloud as compositionModel_(nullptr), and then in the body of this constructor, it is reset by using function setModels(). This is quite similar to the setModels() we used in the constructor body of ReactingMultiphaseCloud, therefore it will not be read through again here. One interesting thing should be noted here is a conversion operator used in autoPtr. As the return type of the function composition() is const Foam::CompositionModel<Foam::ReactingCloud<CloudType>>&, which is not a pointer. We dig into the class autoPtr and find out that it use the conversion operator:

1
2
3
4
5
template<class T>
inline Foam::autoPtr<T>::operator const T&() const
{
    return operator()();
}

When the conversion operator is used, the conversion type should be mentioned after the “operator” sign as something like:

1
operator Type(){}

The conversion operator is used to transform the type of the class to some other types. Read this more in C++ books.

The compositionModel used in the simplifiedSiwek is called singleMixtureFraction. During the reset of the compositionModel_ object, in the related New function, a pointer to the class singleMixtureFraction will be returned. The constructor of singleMixtureFraction is show 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
32
33
34
35
36
template<class CloudType>
Foam::SingleMixtureFraction<CloudType>::SingleMixtureFraction
(
    const dictionary& dict,
    CloudType& owner
)
:
    CompositionModel<CloudType>(dict, owner, typeName),

    idGas_(-1),
    idLiquid_(-1),
    idSolid_(-1),

    YMixture0_(3)
{
    constructIds();

    if (this->phaseProps().size() != 3)
    {
        FatalErrorInFunction
            << "Incorrect number of phases: " << nl
            << "    Please specify 1 gas, 1 liquid and 1 solid"
            << exit(FatalError);
    }

    this->coeffDict().lookup("YGasTot0") >> YMixture0_[idGas_];
    this->coeffDict().lookup("YLiquidTot0") >> YMixture0_[idLiquid_];
    this->coeffDict().lookup("YSolidTot0") >> YMixture0_[idSolid_];

    if (mag(sum(YMixture0_) - 1.0) > small)
    {
        FatalErrorInFunction
            << "Sum of phases should be 1. Phase fractions:" << nl
            << YMixture0_ << exit(FatalError);
    }
}

The function constuctIds() 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
24
25
26
27
28
29
30
31
 template<class CloudType>
 void Foam::SingleMixtureFraction<CloudType>::constructIds()
 {
     forAll(this->phaseProps(), phaseI)
     {
         switch (this->phaseProps()[phaseI].phase())
         {
             case phaseProperties::GAS:
             {
                 idGas_ = phaseI;
                 break;
             }
             case phaseProperties::LIQUID:
             {
                 idLiquid_ = phaseI;
                 break;
             }
             case phaseProperties::SOLID:
             {
                 idSolid_ = phaseI;
                 break;
             }
             default:
             {
                 FatalErrorInFunction
                     << "Unknown phase enumeration" << nl << abort(FatalError);
             }
         }
     }
...
 }

In the constructor of CompositionModel, the phaseProps_ is initialized as:

1
2
3
4
5
6
7
phaseProps_
(
    this->coeffDict().lookup("phases"),
    thermo_.carrier().species(),
    thermo_.liquids().components(),
    thermo_.solids().components()
)

The phaseProps_ is defined as a phasePropertiesList :

1
2
//- List of phase properties
phasePropertiesList phaseProps_;

In the constructor of phasePropertiesList:

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::phasePropertiesList::phasePropertiesList
(
    Istream& is,
    const wordList& gasNames,
    const wordList& liquidNames,
    const wordList& solidNames
)
:
    props_(is),
    phaseTypeNames_(),
    stateLabels_()
{
    forAll(props_, i)
    {
        props_[i].reorder(gasNames, liquidNames, solidNames);
    }

    phaseTypeNames_.setSize(props_.size());
    stateLabels_.setSize(props_.size());
    forAll(props_, i)
    {
        phaseTypeNames_[i] = props_[i].phaseTypeName();
        stateLabels_[i] = props_[i].stateLabel();
    }
}

The props_ is defined as a List object:

1
2
//- List of phase properties
List<phaseProperties> props_;

This List will be initialized by the Istream object 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
44
45
void Foam::phaseProperties::reorder
(
    const wordList& gasNames,
    const wordList& liquidNames,
    const wordList& solidNames
)
{
    // Determine the addressing to map between species listed in the phase
    // with those given in the (main) thermo properties
    switch (phase_)
    {
        case GAS:
        {
            // The list of gaseous species in the mixture may be a sub-set of
            // the gaseous species in the carrier phase
            setCarrierIds(gasNames);
            break;
        }
        case LIQUID:
        {
            // Set the list of liquid species to correspond to the complete list
            // defined in the thermodynamics package.
            reorder(liquidNames);
            // Set the ids of the corresponding species in the carrier phase
            setCarrierIds(gasNames);
            break;
        }
        case SOLID:
        {
            // Set the list of solid species to correspond to the complete list
            // defined in the thermodynamics package.
            reorder(solidNames);
            // Assume there is no correspondence between the solid species and
            // the species in the carrier phase (no sublimation).
            break;
        }
        default:
        {
            FatalErrorInFunction
                << "Invalid phase: " << phaseTypeNames[phase_] << nl
                << "    phase must be gas, liquid or solid" << nl
                << exit(FatalError);
        }
    }
}
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
void Foam::phaseProperties::setCarrierIds
(
    const wordList& carrierNames
)
{
    carrierIds_ = -1;

    forAll(names_, i)
    {
        forAll(carrierNames, j)
        {
            if (carrierNames[j] == names_[i])
            {
                carrierIds_[i] = j;
                break;
            }
        }
        if (carrierIds_[i] == -1)
        {
            FatalErrorInFunction
                << "Could not find carrier specie " << names_[i]
                << " in species list" <<  nl
                << "Available species are: " << nl << carrierNames << nl
                << exit(FatalError);
        }
    }
}
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

void Foam::phaseProperties::reorder(const wordList& specieNames)
{
    // ***HGW Unfortunately in the current implementation it is assumed that
    // if no species are specified the phase is not present and this MUST
    // be checked at the point of use.  This needs a rewrite.
    if (!names_.size())
    {
        return;
    }

    // Store the current sames and mass-fractions
    List<word> names0(names_);
    scalarField Y0(Y_);
    
    // Update the specie names to those given
    names_ = specieNames;
    
    // Re-size mass-fractions if necessary, initialize to 0
    if (names_.size() != names0.size())
    {
        Y_.setSize(names_.size());
        Y_ = 0;
    }
    
    // Set the mass-fraction for each specie in the list to the corresponding
    // value in the original list
    forAll(names0, i)
    {
        bool found = false;
        forAll(names_, j)
        {
            if (names_[j] == names0[i])
            {
                Y_[j] = Y0[i];
                found = true;
                break;
            }
        }
    
        if (!found)
        {
            FatalErrorInFunction
                << "Could not find specie " << names0[i]
                << " in list " <<  names_
                << " for phase " << phaseTypeNames[phase_]
                << exit(FatalError);
        }
    }
}
1