CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Issues in FGM combustion model implementation

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By hk318i

Reply
 
LinkBack Thread Tools Display Modes
Old   September 12, 2015, 16:46
Default Issues in FGM combustion model implementation
  #1
Member
 
Lisandro Maders
Join Date: Feb 2013
Posts: 67
Rep Power: 5
Lisandro Maders is on a distinguished road
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:

Code:
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:
Code:
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:
Code:
 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


Lisandro Maders is offline   Reply With Quote

Old   September 13, 2015, 12:18
Default
  #2
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 219
Rep Power: 10
hk318i is on a distinguished road
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.

Quote:
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.

Quote:
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.

Quote:
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.

Quote:
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 likes this.
hk318i is online now   Reply With Quote

Old   September 13, 2015, 15:35
Default
  #3
Member
 
Lisandro Maders
Join Date: Feb 2013
Posts: 67
Rep Power: 5
Lisandro Maders is on a distinguished road
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 is offline   Reply With Quote

Old   September 18, 2015, 20:51
Default
  #4
Member
 
Lisandro Maders
Join Date: Feb 2013
Posts: 67
Rep Power: 5
Lisandro Maders is on a distinguished road
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:

Code:
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.

Code:
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.

Code:
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
Lisandro Maders is offline   Reply With Quote

Old   September 20, 2015, 06:07
Default
  #5
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 219
Rep Power: 10
hk318i is on a distinguished road
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
hk318i is online now   Reply With Quote

Old   September 21, 2015, 10:33
Default
  #6
Member
 
Lisandro Maders
Join Date: Feb 2013
Posts: 67
Rep Power: 5
Lisandro Maders is on a distinguished road
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:

Code:
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 is offline   Reply With Quote

Old   September 21, 2015, 18:25
Default
  #7
Member
 
Lisandro Maders
Join Date: Feb 2013
Posts: 67
Rep Power: 5
Lisandro Maders is on a distinguished road
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

Code:
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..):


Code:
// * * * * * * * * * * * * * 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
Code:
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:

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:

Code:
       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:

Code:
        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 is offline   Reply With Quote

Old   September 21, 2015, 18:28
Default
  #8
Member
 
Lisandro Maders
Join Date: Feb 2013
Posts: 67
Rep Power: 5
Lisandro Maders is on a distinguished road
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..
Attached Images
File Type: png Screenshot from 2015-09-21 19:13:28.png (38.9 KB, 13 views)
File Type: png Screenshot from 2015-09-21 19:14:17.png (12.7 KB, 10 views)
Lisandro Maders is offline   Reply With Quote

Old   September 29, 2015, 10:30
Default
  #9
Member
 
Lisandro Maders
Join Date: Feb 2013
Posts: 67
Rep Power: 5
Lisandro Maders is on a distinguished road
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
Lisandro Maders is offline   Reply With Quote

Reply

Tags
constructor, fgm, inheritance

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Diff in Combustion Modelling By Species Transport Model & Non-Premixed Combustion bhanuday.sharma FLUENT 0 June 23, 2015 02:50
Liquid rocket engine non-premixed combustion model Erik FLUENT 2 November 28, 2013 15:09
eddy dissipation model: combustion doesn't occur roukaia FLUENT 0 December 24, 2011 10:10
the combustion model problem hanshandiaoke CD-adapco 0 November 25, 2011 21:37
combustion model Hennie van der Westhuizen CD-adapco 7 February 27, 2002 03:10


All times are GMT -4. The time now is 21:56.