[Lisandro Maders]

Hello,

I will divide this text a bit in order to organize the information:

1. If you want to be situated in the project I am doing
You can find information about the FGM combustion model at: http://www.fgm-combustion.org/
My final project consists in implement such method to solve a non-premixed laminar jet flame (2D domain). The deadline is 1st of November and I am struggling a bit to understand some parts of a code I have.

2. Basic idea of implementation of FGM in OpenFOAM
The base idea is to solve 2 additional transport equations instead of specie transport and energy equation (in my case, but there are other dimensions of manifold to solve), and update the diffusion coefficients reading a table (called manifold). The properties I need to read are:
- viscosity for N-S equation
- density for continuity and N-S equation
- Cp and thermal conductivity for the 2 additional transport equations

Such values are parametrized with a combination of the other 2 added transport variables and read from table according to the such transport variables values at table.
This table is created by 1D simulations varying the stretch rate flame value to take int account different curvatures of the flame in the 2D simulation.

3. Problems in OpenFOAM understanding
I got some files of a FGM implementation done in OpenFOAM and I am trying to debug some parts of the code. I will start putting the difficulties I am facing to here.

3.1 basicPsiThermo class in OF 2.3.0
It uses as base class the basicPsiThermo class. I am using OpenFOAM 2.3.0 version and there is not such class on it. I read in a thread that it's name have been changed to psiThermo however that was not an official statement. I could't find in OpenFOAM releases anything talking about such change.

3.2 Inheritance concept
Analysing the following pieces of code..

hPsiFanzy.H file:

class hPsiFanzy
:
    public basicPsiThermo
{
    // Private data

        //- Sensible enthalpy field [J/kg]
        volScalarField h_;

        //- Table look up routines
        const fanzyLookUp& fgmTable_;
        
        //- Reference to CV1 & CV2
        const volScalarField& foamCV1_;
        const volScalarField& foamCV2_;
        
        //- FGM table index maps (CpIndex, RIndex, muIndex, alphaIndex)
        const labelList fgmThermoTransportIndices_;
        
        
    // Private Member Functions

        //- Calculate the thermo variables
        void calculate();

        //- Construct as copy (not implemented)
        hPsiFanzy(const hPsiFanzy&);


public:

    //- Runtime type information
    TypeName("hPsiFanzy");

    // Constructors

        //- Construct from mesh
        hPsiFanzy
        (
            const fvMesh&,
            const fanzyLookUp&,
            const volScalarField&,
            const volScalarField&
        );


    //- Destructor
    ~hPsiFanzy();

etc..

Ok, here the class hPsiFanzy is declared as being a derived class of basicPsiThermo class, and have h_ , fgmTable_ , foamCV1_ , foamCV2_ , fgmThermoTransportIndices_ as private members and calculate() as private member function..

Then in the hPsiFanzy.C file, the definition of the constructor is as follows:

Foam::hPsiFanzy::hPsiFanzy
(
    const fvMesh& mesh,
    const fanzyLookUp& fanzyLookUp,
    const volScalarField& foamCV1,
    const volScalarField& foamCV2
)
:
    basicPsiThermo(mesh),

    h_
    (
        IOobject
        (
            "h",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimEnergy/dimMass,
        this->hBoundaryTypes()
    ),

    fgmTable_(fanzyLookUp),
    foamCV1_(foamCV1),
    foamCV2_(foamCV2),
    
    fgmThermoTransportIndices_(lookup("fgmThermoTransportIndices"))
{
.... //implementation of functionalities of constructor

My big concern is that I thought that in the constructor definition of a derived class it should be passed as parameter only constructors of the base class after the : sign.
That's what I meant it should be a standard derived class definition:

 class base
          {
         protected:
                   int i;
         public:
                base(int x) { i = x; }
                ~base() {}
         };
   class derived : public base
          {
               int j;
         public:
                derived(int x, int y) : base(y) { j = x;}
                ~derived() {}
         };

But as you can see in red the private members from the derived class are being passed. I understand it should be doing this in order to the base class can use such members, but I thought it should be written by constructing a basiPsiThermo object with such arguments, not only putting the arguments like this.

Another thing that bothers me is that in the hPsiFanzy.H file the object h_ is declared as being from volScalarField class and has been constructed with no arguments (see blue). Then later in the constructor definition it has been constructed with IOobject argument. I was not able to understand that, maybe I am misunderstanding a lot some concepts..

Briefing, I would like to understand the concept behind this hPsiFanzy as a class, because it is one of the most important classes in the FGM implementation. Understanding it I will be able to start understanding the whole implementation..

3.3 variable_ vs. variable
Is there some meaning of naming a variable like this variable_ instead of variable in OpenFOAM? I mean, in general it is used to identify some type of variable or usability of it?

Thanks for any reply! I really need such help with this.


Best regards,

Lisandro

[hk318i]

Hello,

Most of your questions are more related to C++ programming concepts more than they are related to OpenFOAM.
I will try to give you concise answers and clues to your questions.

3.1 basicPsiThermo class in OF 2.3.0
It uses as base class the basicPsiThermo class. I am using OpenFOAM 2.3.0 version and there is not such class on it. I read in a thread that it's name have been changed to psiThermo however that was not an official statement. I could't find in OpenFOAM releases anything talking about such change.

That is true there is no class called basicPsiThermo starting from version-2.2, and I don't know exactly what are the changes. You could compare thermophysicalModels library from version-2.3 against version-2.2. Or even better compare thermophysicalModels library from version-2.3 against the version which your code depends on originally. It should be straight forward to construct more or less the same basicPsiThermo model using the available template classes in OpenFOAM-2.3.

3.2 Inheritance concept
Analysing the following pieces of code..

My big concern is that I thought that in the constructor definition of a derived class it should be passed as parameter only constructors of the base class after the : sign.
That's what I meant it should be a standard derived class definition:

Not true, All the private data members should be initialised after initialising the base class.

3.2 Inheritance concept
Another thing that bothers me is that in the hPsiFanzy.H file the object h_ is declared as being from volScalarField class and has been constructed with no arguments (see blue). Then later in the constructor definition it has been constructed with IOobject argument. I was not able to understand that, maybe I am misunderstanding a lot some concepts..

It is a common practice to declare the data and functions in the header (.H) file and initialise the data in the constructor and define the functions in the source file (.C). About h_ in particular, volScalarField could be constructed from IOobject.

3.3 variable_ vs. variable
Is there some meaning of naming a variable like this variable_ instead of variable in OpenFOAM? I mean, in general it is used to identify some type of variable or usability of it?

Generally speaking there is no difference BUT it is a practice which makes the code easier to follow and maintain. variable_ means a data member in the class which is accessible by any function in the class. variable is a local variable within the limited scope of a member function. It is OpenFOAM style, however different projects and developers have their own style. But generally, they are always trying to differentiate between class data members and local variables.

I hope that, I gave you a hint about your problem and good luck with your project.

Best wishes,
Hassan

[Lisandro Maders]

Dear Hassan,

thanks a lot! I was able to understand everything and also check some c++ concepts in order to ensure I have understood. Thank you very much, your answer was clear and good.

Best regards,

Lisandro

[Lisandro Maders]

Hi again,

now I am struggling a bit to understand the logic of OpenFOAM fields updating.

What I want to do: I want the following properties to be read from a table --> rho, mu, alpha and Cp. Therefore, I do not need to calculate such properties using a transport model, equation of state or thermodynamics model as it is in most cases.

The problem: let's say that the way I started doing it was creating a thermophysical model based in some existent one (like psiThermo or rhoThermo) and did it by deriving my model from one of these. Let's say I am using the psiThermo as example and I want to debug how the Cp is being calculated.

The protected members of psiThermo are written below, at psiThermo.H file:

class psiThermo
:
    public fluidThermo
{

protected:

    // Protected data

        //- Compressibility [s^2/m^2]
        volScalarField psi_;

        //- Dynamic viscosity [kg/m/s]
        volScalarField mu_;

As you can see, there is no member for Cp, like Cp_. Also, there is no member function calling Cp. So, let's look at fluidThermo, which is the class from which psiThermo is derived.

class fluidThermo
:
    public basicThermo
{

Again, there is no private or protected member called Cp neither member function Cp(). So, let's take a look at basicThermo, which is the base class of fluidThermo.

class basicThermo
:
    public IOdictionary
{

protected:

    // Protected data

        //- Phase-name
        const word& phaseName_;


        // Fields

            //- Pressure [Pa]
            volScalarField& p_;

            //- Temperature [K]
            volScalarField T_;

            //- Laminar thermal diffusuvity [kg/m/s]
            volScalarField alpha_;

......
.......
......
......

            //- Heat capacity at constant pressure [J/kg/K]
            virtual tmp<volScalarField> Cp() const = 0;

            //- Heat capacity at constant pressure for patch [J/kg/K]
            virtual tmp<scalarField> Cp
            (
                const scalarField& p,
                const scalarField& T,
                const label patchi
            ) const = 0;
...
....
....
.....

Ok, now I can find some Cp in the code. However, it is not a member, it is a funcion that returns a volScalarFied, what I would suppose it is the Cp_ field. However, as you can see it only initializes it to 0. Then you cannot find Cp_ anywhere else in the code.

So my question is, I can imagine that the Cp is calculated through the thermodynamic model you choose at the thermoPhysicalProperties dictionary like janaf or hConst, but where is it being called?

And thinking about my problem, that I do not need to calculate such values (only need to read them at a table and use them in the diffusion terms of transport equations) which is the best way to implement it?


Thank you in advance for any attempt of help!


Regards,

Lisandro

[hk318i]

Hello,

The thermodynamic library in OF is really hard to follow. It has a very clever implementation but it is not straightforward. Unfortunately, I cannot give you a definitive answer to your question. You conclusion is quiet right but please note that hConst and janaf are template classes, therefore the cp values are calculated not only based on hConst or janaf but also depend on equationOfState.

best wishes,
Hassan

[Lisandro Maders]

Yes I understand it is not that trivial.. Even if what I post above is quiet right, it is still not enough for me to understand properly how can I deal with such fields in order to make them to be read from a table and use such updated values in equations instead of using the thermodynamics and transport models (and equation of state for rho)..

What I have done so far is to right a new thermophysical model based on basicPsiThermo, and tried to explicitly assign the table values to the properties trough commands like this:

rhoCells[celli] = myTable_.getValue2D
                         (
                             inputVar1Cells[celli],
                             inputVar2Cells[celli],
                             fgmThermoTransportIndices_[6]
                         );

Where the number 6 is the correspondent index in Table that contains rho.

Of course I have instantiated a rho_ object since the basicPsiThermo did not have it..

What the function getValue2D does is only to read 2 other variables in table, make a combination of their values, interpolate and then take the correspondent rho value in table. The function is 100% sure working fine (since it was working for other properties in another solver), but when I plot the rho in paraView I cannot see any variations of my initial condition, what makes me think that it is not properly reading or writing such field. Also, my big doubt is whether the equations are using the rho from table or still using the equationOfState calculation.

I know my post is very general, but if anyone has experience in dealing with fields or, even better, have already experienced working with tables for properties reading, and are available to give some tips to me I can clarify the post.

Thanks a lot,

Lisandro

[Lisandro Maders]

If someone wants to help me a lot so just try to answer this following question, then I will be able to figure out the other stuff..

I have a thermophysical model class called hPsiNew. It is derived from basicPsiThermo (it is from OF2.1.0), which is derived from basicThermo class. Remembering that what I want is to read some fields from a table, I have the following code for hPsiNew class:

1. hPsiNew.H file

class hPsiNew
:
    public basicPsiThermo
{
        const fanzyLookUp& fgmTable_;
        
        const volScalarField& foamCV1_;
        const volScalarField& foamCV2_;
        
        const labelList fgmThermoTransportIndices_;
                
    // Private Member Functions

        void calculate();

        //- Construct as copy (not implemented)
        hPsiNew(const hPsiNew&);

public:

    TypeName("hPsiNew);

        //- Construct from mesh
        hPsiNew
        (
            const fvMesh&,
            const fanzyLookUp&,
            const volScalarField&,
            const volScalarField&
        );


    //- Destructor
    ~hPsiFanzy();

        void correct();

            tmp<scalarField> Cp
            (
                const scalarField& T,
                const label patchi
            ) const;

            tmp<volScalarField> Cp() const;
};

}

and my hPsiNew.C file (I will put only what matters..):

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

void Foam::hPsiFanzy::calculate()
{
    scalarField& TCells = this->T_.internalField();
    scalarField& psiCells = this->psi_.internalField();
    scalarField& muCells = this->mu_.internalField();
    scalarField& alphaCells = this->alpha_.internalField();

    const scalarField& foamCV1Cells = foamCV1_.internalField();
    const scalarField& foamCV2Cells = foamCV2_.internalField();
    
    forAll(TCells, celli)
    {
        scalar Ttab = fgmTable_.getValue2D
                      (
                          foamCV1Cells[celli],
                          foamCV2Cells[celli],
                          fgmThermoTransportIndices_[4]
                      );
        TCells[celli] = Ttab;
        
        // psi = 1.0/(R*T);
       psiCells[celli] = 1.0/(TCells[celli]*fgmTable_.getValue2D
                                            (
                                                foamCV1Cells[celli],
                                                foamCV2Cells[celli],
                                                fgmThermoTransportIndices_[1]
                                            ));
    }
    forAll(T_.boundaryField(), patchi)
    {
        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];

        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];

        const fvPatchScalarField& pfoamCV1 = foamCV1_.boundaryField()[patchi];
        const fvPatchScalarField& pfoamCV2 = foamCV2_.boundaryField()[patchi];
        
        if (pT.fixesValue())
        {
            forAll(pT, facei)
            {
                scalar Ttab = fgmTable_.getValue2D
                                (
                                    pfoamCV1[facei],
                                    pfoamCV2[facei],
                                    fgmThermoTransportIndices_[4]
                                );

               ppsi[facei] = 1.0/(pT[facei]*fgmTable_.getValue2D
                                               (
                                                   pfoamCV1[facei],
                                                   pfoamCV2[facei],
                                                   fgmThermoTransportIndices_[1]
                                               ));
            }
        }



...
...
...

// Constructor

Foam::hPsiFanzy::hPsiFanzy
(
    const fvMesh& mesh,
    const fanzyLookUp& fanzyLookUp,
    const volScalarField& foamCV1,
    const volScalarField& foamCV2
)
:
    basicPsiThermo(mesh),

    fgmTable_(fanzyLookUp),
    foamCV1_(foamCV1),
    foamCV2_(foamCV2),
    
    fgmThermoTransportIndices_(lookup("fgmThermoTransportIndices"))
{
    scalarField& TCells = this->T_.internalField();

    const scalarField& foamCV1Cells = foamCV1_.internalField();
    const scalarField& foamCV2Cells = foamCV2_.internalField();

    forAll(TCells, celli)
    {
        // Initialize enthalpy and temperature with values from the FGM table
        TCells[celli] = fgmTable_.getValue2D
                        (
                            foamCV1Cells[celli],
                            foamCV2Cells[celli],
                            fgmThermoTransportIndices_[4]
                        );
    }

...
...
...

    forAll(T_.boundaryField(), patchi)
    {
        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
        
        const fvPatchScalarField& pfoamCV1 = foamCV1_.boundaryField()[patchi];
        const fvPatchScalarField& pfoamCV2 = foamCV2_.boundaryField()[patchi];
            
        // Initialize enthalpy and temperature with values from the FGM table
        if (pT.fixesValue())
        {
            forAll(pT, facei)
            {
                scalar Ttab = fgmTable_.getValue2D
                                (
                                    pfoamCV1[facei],
                                    pfoamCV2[facei],
                                    fgmThermoTransportIndices_[4]
                                );
            }
        }
        else
        {
            forAll(pT, facei)
            {
                pT[facei] = fgmTable_.getValue2D
                                (
                                    pfoamCV1[facei],
                                    pfoamCV2[facei],
                                    fgmThermoTransportIndices_[4]
                                );
            }
        }
    }

...
...
...

So my questions are:

1. T_ object
The first line of calculate() member function does this

scalarField& TCells = this->T_.internalField();

The thing is, this T_ object is only declared and defined at basicThermo class. That's ok, it is using such class to deal with T fields.. Then, in the first loop we have these lines of code:

scalar Ttab = fgmTable_.getValue2D
                      (
                          foamCV1Cells[celli],
                          foamCV2Cells[celli],
                          fgmThermoTransportIndices_[4]
                      );
        TCells[celli] = Ttab;

So I assume that from now on the TCells (the T value of the Cells) are going to be always the T values from the table. So it will never more use the basicThermo class to deal with T field, right?

2. psi field
As you can see, the psi field is calculated by 1/(R*T), being both R and T values taken from table.
The thing is:
a) this is the psi field when I solved the hEqn together with the other equations. It means it was calculating some T field (but I am not sure which one was used to calculate psi). These area the results for psi (at some time -not important)

See FIgure a) in next post

b) Then, when I only deactivate the hEqn these are the results for psi (at the same time):

See Figure b) in next post

So, my conclusion is that it is not properly doing the following command:

       psiCells[celli] = 1.0/(TCells[celli]*fgmTable_.getValue2D
                                            (
                                                foamCV1Cells[celli],
                                                foamCV2Cells[celli],
                                                fgmThermoTransportIndices_[1]
                                            ));

The Tcells[celli] should be the T read from table, as it is told to do at:

        scalar Ttab = fgmTable_.getValue2D
                      (
                          foamCV1Cells[celli],
                          foamCV2Cells[celli],
                          fgmThermoTransportIndices_[4]
                      );
        TCells[celli] = Ttab;

But apparently it is not, since when I deactivate the hEqn the psi field changed..

So I am struggling to deal with field's writting and reading stuff... If anyone could help me on this it would be great!


Best regards,

Lisandro

[Lisandro Maders]

Figure a) is the most non-uniform one

Figure b) is the most uniform one


Obs: I struggled to put images I do not know why,sorry about this confusion..

[Lisandro Maders]

I solved my problem regarding to properties being taken from table. It is working fine now. If anyone have troubles with this in future ask me then I can give a help. I will not post here because I was doing so many things wrong but now I understood how to do this and a bit more how to deal with fields in OpenFOAM.

Lisandro

[R.Mah]

Hallo Lisandro Maders
Now I am working on the implementation of the FGM combustion model into OpenFOAM and I would really appreciate if you can help on that by posting the code that you worked on or if you know another code with more recent OF version.
Thanks

[spk363]

I would also like to know if there are implementations of the FGM model in newer versions like OF 4.1.

[Lisandro Maders]

Originally Posted by R.Mah

Hallo Lisandro Maders
Now I am working on the implementation of the FGM combustion model into OpenFOAM and I would really appreciate if you can help on that by posting the code that you worked on or if you know another code with more recent OF version.
Thanks

Hi R.Mah,

I coded the FGM model for non-premixed laminar flames in OpenFOAM 2.1.x and since then I have been working with this version when I want to use FGM. I did not put effort to adapt it to newer versions.

If you want I can share the code with you. It was the first big code I did in OpenFOAM so it is not perfectly programmed in terms of OpenFOAM style.

Let me know if such solver for non-premixed laminar flames suits for you.

Regards, Lisandro Maders

[Lisandro Maders]

Originally Posted by Lisandro Maders

Hi R.Mah,

I coded the FGM model for non-premixed laminar flames in OpenFOAM 2.1.x and since then I have been working with this version when I want to use FGM. I did not put effort to adapt it to newer versions.

If you want I can share the code with you. It was the first big code I did in OpenFOAM so it is not perfectly programmed in terms of OpenFOAM style.

Let me know if such solver for non-premixed laminar flames suits for you.

Regards, Lisandro Maders

One important thing to state is that this solver considers that you have already created a manifold (the lookup table). So it only solves the multidimensional step of FGM solution. The solution of the flamelets (1D flames) and the creation of the manifold are not accounted for in this solver. Our group (Combustion Laboratory at UFRGS) uses a in-house code to solve the flamelets and create the manifolds.

Regards,

Lisandro Maders

[Tobi]

Dear Lisandro,


I re-activate this thread as I want to know a few things here which you might be able to answer:

• Does the FGM Solver exist on a repository? If yes, would you be so kind to share it or is it (as mostly) an inhouse-code
• The flamelet generator is an inhouse code. Would you be so kind to tell me why everybody is doing is thermo-kinetic calculation with closed codes? I only know Cantera which can perform kinetic calculations; however, for flamelet models we might solve directly for the mixture fraction rather than a space distance
• Can you give an insight into the differences between the laminar flamelet model introduced by Peters and the flamelet generalized manifold approach? I am not familiar with the latter one.

Furthermore, do you know anybody who would like to create an open-source kinetic solver for, e.g., flamelets?

Kind regards and thanks in advance,
Tobias

[zhangyan]

Originally Posted by Tobi

Dear Lisandro,


I re-activate this thread as I want to know a few things here which you might be able to answer:

• Does the FGM Solver exist on a repository? If yes, would you be so kind to share it or is it (as mostly) an inhouse-code
• The flamelet generator is an inhouse code. Would you be so kind to tell me why everybody is doing is thermo-kinetic calculation with closed codes? I only know Cantera which can perform kinetic calculations; however, for flamelet models we might solve directly for the mixture fraction rather than a space distance
• Can you give an insight into the differences between the laminar flamelet model introduced by Peters and the flamelet generalized manifold approach? I am not familiar with the latter one.

Furthermore, do you know anybody who would like to create an open-source kinetic solver for, e.g., flamelets?

Kind regards and thanks in advance,
Tobias

Dear Tobias,
There is a widely-used open source 1D chemical calculator, called FlameMaster.
You can find it here: https://www.itv.rwth-aachen.de/index.php?id=13
FlameMaster can solve 1D flamelets equations on the mixture fraction coordinate.
And another one from a personal developer:
https://github.com/ZX114/ChemFlow

And some related repositories:
https://github.com/weixian001/FPVFoam-v3.0.1
https://github.com/CCChuanfeng/premixedFGMFoam
https://github.com/ZX114/sprayFlameletFoam

[blttkgl]

Just to chime in, Ember is also another open-source unsteady flamelet solver, developed by the same people developing Cantera. I have been testing its performance with respect to Chem1D (another closed-source flamelet solver, however possible to get from TU Eindhoven), and it looks quite promising.


https://github.com/speth/ember


Bulut