[kaizhangqmul]

Hi,

I am trying to output fuel reaction rate defined on each cell for further calculations (reactingFoam). What I did is :
1) In createFields.H, add :
autoPtr<psiChemistryModel> pChemistry(psiChemistryModel::New(mesh));
psiChemistryModel& chemistry = pChemistry();
volScalarField::Internal rr
(
IOobject
(
"rr",
runTime.timeName(),
mesh,
IOobject::NO_READ, //READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh
);
2) In myreactingFoam.C, add:
rr = chemistry.RR(2);

Compilation has no problem, but as you can see in log attached:
reaction rate = dimensions [1 -3 -1 0 0 0 0]
value uniform 0; --->I am expecting the values calculated from Arrenius law based on cell mass fraction which shouldn't be 0.

Can anyone help pls?

Thanks in advance.
Kai.

[adhiraj]

I was able to access the reaction rates in OpenFOAM 2.3.x like this:

    forAll(Y, i)
    {
        ROP_[i].field() = -(reaction->R(Y[i])()).source()/mesh.V();
    }

Maybe you can see if there is an equivalent for version 5.x?

[kaizhangqmul]

Hi,

Thanks for your reply. I will see if I could find similar functions. Btw, did you add anything in createFields.H such as how is 'reaction' defined?( something like what I have done above?)

Kai.

[adhiraj]

Unfortunately at the moment I am not familiar enough with the source code of OpenFOAM 5.x to tell you the exact changes you need to make.

In version 2.3.x for reactingFoam, the variable reaction is defined at the top of the createFields.H. This is already there as part of OpenFOAM's original code:

Info<< "Creating reaction model\n" << endl;

autoPtr<combustionModels::psiCombustionModel> reaction
(
    combustionModels::psiCombustionModel::New(mesh)
);

I added a variable to store and print the reaction rates like this:

PtrList<volScalarField> ROP_(Y.size());
forAll(Y, i)
{
    const word RRname = "ROP." + Y[i].name();
    ROP_.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                RRname,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
        )
    );    
}

[kaizhangqmul]

Hi,

Many thanks, it seems that also works in 5.x. I will to some more checks to make sure it does return the right values.

Thanks again!
Kai.

[faeze.d]

Originally Posted by adhiraj

I was able to access the reaction rates in OpenFOAM 2.3.x like this:

    forAll(Y, i)
    {
        ROP_[i].field() = -(reaction->R(Y[i])()).source()/mesh.V();
    }

Maybe you can see if there is an equivalent for version 5.x?

Hi, Thanks for the post. It was very helpful.
I think we should not divide R by mesh.V
Coz the dimension of RR in chemistryModel.C is kg/(m3s). Don't ya think so?

[zhangyan]

Originally Posted by faeze.d

Hi, Thanks for the post. It was very helpful.
I think we should not divide R by mesh.V
Coz the dimension of RR in chemistryModel.C is kg/(m3s). Don't ya think so?

Yes, RR in chemistryModel.C is kg/(m3s).
But, R(volScalarField& Y) in combustionModel class returns an fvScalarMatrix, which is kg/s.
Ref:Chemistry Source Term RR()

[kaizhangqmul]

Originally Posted by faeze.d

Hi, Thanks for the post. It was very helpful.
I think we should not divide R by mesh.V
Coz the dimension of RR in chemistryModel.C is kg/(m3s). Don't ya think so?

Yes, you are right regarding the RR unit = kg(m3s), but when we want to use it in equation or further calculation, what we need will be volumetric RR, i.e RR divided by mesh.V.

[faeze.d]

Originally Posted by kaizhangqmul

Yes, you are right regarding the RR unit = kg(m3s), but when we want to use it in equation or further calculation, what we need will be volumetric RR, i.e RR divided by mesh.V.

Thanks for your prompt reply. I'm confused because in Yi equation we have:
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->muEff(), Yi)
==
reaction->R(Yi)
+ fvOptions(rho, Yi)
the dimension of all terms are kg/(m3s). so the source term should have the same dimension.
I don't have a clue why we divide it by volume again?

[kaizhangqmul]

Originally Posted by faeze.d

Thanks for your prompt reply. I'm confused because in Yi equation we have:
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
- fvm::laplacian(turbulence->muEff(), Yi)
==
reaction->R(Yi)
+ fvOptions(rho, Yi)
the dimension of all terms are kg/(m3s). so the source term should have the same dimension.
I don't have a clue why we divide it by volume again?

My mistake, I was thinking of an energy equation in my code, and get confused with the code here.
The reaction->R() below has a unit of kg/s, and what you need will be R/mesh.V in your Yi equation.

forAll(Y, i)
{
ROP_[i].field() = -(reaction->R(Y[i])()).source()/mesh.V();
}

[felixFoam]

Originally Posted by adhiraj

Unfortunately at the moment I am not familiar enough with the source code of OpenFOAM 5.x to tell you the exact changes you need to make.

In version 2.3.x for reactingFoam, the variable reaction is defined at the top of the createFields.H. This is already there as part of OpenFOAM's original code:

Info<< "Creating reaction model\n" << endl;

autoPtr<combustionModels::psiCombustionModel> reaction
(
    combustionModels::psiCombustionModel::New(mesh)
);

I added a variable to store and print the reaction rates like this:

PtrList<volScalarField> ROP_(Y.size());
forAll(Y, i)
{
    const word RRname = "ROP." + Y[i].name();
    ROP_.set
    (
        i,
        new volScalarField
        (
            IOobject
            (
                RRname,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
        )
    );    
}

Hello adhiraj,

I tried your solution, added your code to createFields.H and compiled the solver with wmake. However, if I run a case with my solver, the simulation works but the reaction rates are not printed out. Did I miss something?

Regards

[adhiraj]

Originally Posted by felixFoam

Hello adhiraj,

I tried your solution, added your code to createFields.H and compiled the solver with wmake. However, if I run a case with my solver, the simulation works but the reaction rates are not printed out. Did I miss something?

Regards

Can you give some more information?

Which version of OpenFOAM is this? Are the variables not written to the time directories at all, or do you get zeros?

[felixFoam]

Originally Posted by adhiraj

Can you give some more information?

Which version of OpenFOAM is this? Are the variables not written to the time directories at all, or do you get zeros?

Thanks a lot for the quick reply,

I was able to fix my issue, it was a small bug somewhere else that I just found by chance. However, thank you a lot, this forum and especially you were incredibly helpful!

Anyway, as this forum left one minor uncertainty (at least for beginners in C++), this might help some future readers:

1. Definition of new field(s) in createFields.H
2. Calculation of the source terms in YEqn.H

(It might work differently, but this is how I implemented it and it sort of makes sense)