CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   Species Reaction Rate RR() (https://www.cfd-online.com/Forums/openfoam/203367-species-reaction-rate-rr.html)

kaizhangqmul June 26, 2018 00:50

Species Reaction Rate RR()
 
1 Attachment(s)
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 June 26, 2018 11:30

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

Code:

    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 June 26, 2018 12:12

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 June 26, 2018 13:13

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:

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:

Code:

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 June 26, 2018 22:02

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 October 16, 2018 03:17

Quote:

Originally Posted by adhiraj (Post 697323)
I was able to access the reaction rates in OpenFOAM 2.3.x like this:

Code:

    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 October 16, 2018 04:57

Quote:

Originally Posted by faeze.d (Post 710143)
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:https://www.cfd-online.com/Forums/op...e-term-rr.html

kaizhangqmul October 16, 2018 05:02

Quote:

Originally Posted by faeze.d (Post 710143)
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 October 16, 2018 05:10

Quote:

Originally Posted by kaizhangqmul (Post 710158)
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 October 16, 2018 05:23

Quote:

Originally Posted by faeze.d (Post 710159)
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 September 12, 2021 12:09

No output
 
Quote:

Originally Posted by adhiraj (Post 697339)
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:

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:

Code:

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 September 20, 2021 15:54

Quote:

Originally Posted by felixFoam (Post 812039)
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 September 21, 2021 06:16

Thank for the quick reply
 
Quote:

Originally Posted by adhiraj (Post 812616)
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)


All times are GMT -4. The time now is 05:44.