[Mirza8]

Hello everyone,

I am running a combustion simulation by the reactingFoam solver in openFOAM 6.0, and everything works fine. I am using a detailed reaction mechanism with 62 species and 262 reactions. After the calculations are finished, I want to write the output and observe the fields of the two following parameters:

1- (My major problem) The rate of each reaction in my domain (262 different fields).
2- (Minor problem) The source term of each species in my domain (62 different fields).

For the first parameter, when I try to look into code to find out where the Arrhenius reaction rates are calculated, I get confused between all the reactionModels, chemistryModels and ODESolvers. Hence I do not know where the rate of each reaction is calculated and how I can get access to it in the solver code. Is there any postProcess tool for this purpose? Or is there any simple code which can be added to the solver to access and get output from the desired fields? I searched the google and the forums but I could not find anything helpful.

However, for the second parameter, we have access to source term of species Y[i] through

reaction->R(Y[i])

which is of type fvMatrix and it is not possible to directly get an output from it. Therefore, someone in THIS TOPIC suggested that we can convert the fvMatrix to volScalarField by multiplying to the respective Y[i] field so something like this:

Ri_[i] = reaction->R(Y[i]) & Y[i];

Where Ri_ is a PtrList<volScalarField>, and it is possible to write out the Ri_[i] fields. I tried this solution and it gives me no error. But can someone verify this to be the correct method and explain why we have to multiply the R(Y[i]) to Y[i]? Isn't this changing the values of the source term fields?

I highly appreciate any suggestion or advice on this matter!

Best regards,
Morteza

[Paulo SPS]

Hi there,
I have the same problem.
I'm still a beginner in openfoam (and C++ for that matter) but I'm not sure the line of code you have there is multiplication.
I think it is bitwise AND between a pointer to reaction.R(Yi) and Y(i), not sure that helps.


I think what you need is something like it is done in the createFields.H for Qdot (I haven't tried it yet):

volScalarField Qdot
(
    IOobject
    (
        "Qdot",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
);

but for each species.


No idea how to get the rate of each reaction though and I'm very interested if you find out how

[Mirza8]

Hi Paulo,

Thanks for your response.

I want to apply your suggestion, but still I do not get how to convert fvMatrix to volScalarField? You see, I did something similar to the code you provided from createFields.h , to create the Ri_ fields:

PtrList<volScalarField> Ri_;
Ri_.setSize(Y.size());

forAll(Y, i)
{
	Ri_.set
	(
		i,
		new volScalarField
		(
			IOobject
			(
				"R_" + Y[i].name(),
				runTime.timeName(),
				mesh,
				IOobject::NO_READ,
				IOobject::AUTO_WRITE
			),
			mesh,
			dimensionedScalar("Ri", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0)
			//reaction->R(Y[i]) & Y[i]
		)
	);
	// if I write Ri_[i] = reaction->R(Y[i]); the compiler complains
	Ri_[i] = reaction->R(Y[i]) & Y[i]; 
	Ri_[i].write(); 
}

But the problem is how you want to convert reaction->R(Y[i]) which is a fvMatrix to volScalar field?? For this conversion, I used the "Ri_[i] = reaction->R(Y[i]) & Y[i];", but I do not know what exactly it does. I think you are right that & is a bitwise AND operator in C++ , but it seems that it also serves as "inner product operator" in OpenFOAM. Just take a look at the attached image which is the operators table from p22 of OpenFOAM's Programmers guide.

So, in case we agree on the fact that & is the inner product operator here, what does "reaction->R(Y[i]) & Y[i];" mean exactly? Is this the correct way to get the source terms for species?

[Paulo SPS]

Hi there,
I didn't know the "&" operator could be used as inner product in openfoam (good to know actually ).


I looked at the code suggested in the topic you refered (Qdot function in singleStepCombustion.C) and it looks like this:

singleStepCombustion<ReactionThermo, ThermoType>::Qdot() const
{
    const label fuelI = singleMixturePtr_->fuelIndex();
    volScalarField& YFuel =
        const_cast<volScalarField&>(this->thermo().composition().Y(fuelI));

    return -singleMixturePtr_->qFuel()*(R(YFuel) & YFuel);
}

and (R(YFuel) & YFuel) does seem to return the rate, but I'll be honest and say I don't understand why (not sure how the structure of fvMatrix works here for the product).

Maybe the best way to be sure is to set a simple case in chemFoam and compare if another software where you can access the rates (eg Cantera).


I also found this topic that handled the problem differently, haven't tried the solution yet though.

[Mirza8]

Thanks again Paul! I have added the other method from the topic that you mentioned to calculate the species source terms. So, now I have two volScalarField variables to calculate the source terms with the following two methods:

// method 1
Ri_[i] = reaction->R(Y[i]) & Y[i];
// method 2
Ri2_[i].field() = -(reaction->R(Y[i])()).source()/mesh.V();

Surprisingly enough, both methods give the same results, with only minor differences at some regions!! For instance, I have plotted the O2 source term calculated by both methods, over a certain line in one of my cases. As you can see in the attached image, the two curves perfectly match each other except in the area close to the inlet (z=0)!

Now I think I can trust the results that we obtain by either of the two mentioned methods! However, it would be nice if somebody could explain why we multiplied reaction->R(Y[i]) to Y[i] in the first method? And why is there a negative sign before (reaction->R(Y[i])()).source() in the second method?

Finally, back to my major problem, can somebody explain how we can get the reaction rates instead of species source terms?

[Mirza8]

I think I have found a solution for my first problem!!!

I have created a modified createFields.H just for post-processing purpose, and I added the following code to it, in order to identify the contribution of each reaction to production of a specific species. For each reaction , this code creates a scalar field which contains the production rate of species through that reaction:

autoPtr<BasicChemistryModel<psiReactionThermo>> tempChem
(
    BasicChemistryModel<psiReactionThermo>::New(thermo)
);

// label (index) of the species of interest
label si = ind_SPEC;

for (label ri=0; ri<tempChem->nReaction(); ri++)
{
	// Calculate the source term for species si, from reaction ri [kg/m^3/s]
	volScalarField::Internal RR
	(
		tempChem->calculateRR(ri, si)
	);
	
	volScalarField RRfield_
	(
		IOobject
		(
			"RR_" + std::to_string(ri),
			runTime.timeName(),
			mesh,
			IOobject::NO_READ,
			IOobject::AUTO_WRITE
		),
		mesh,
		dimensionedScalar("RR", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0)
	);
	RRfield_.field() = RR;
	RRfield_.write();
}

I believe this code gives the thermochemical reaction rates, hence the effects of turbulent combustion models (such as EDC) on reaction rates are not considered. Therefore, the summation of RR(ri, si) over all reactions for a specific species, , is not equal to Ri that we have calculated before by "Ri_[i] = reaction->R(Y[i]) & Y[i];"! Because the later is derived from combustion model and is usually lower than chemical reaction rate!

The code runs and gives the results which qualitatively look fine! However, since I have no data to verify the quantities, I would appreciate it if someone could explain whether I am using the right code or not!

I hope it helps the others who had the same problem as me!

[run]

Hi Mirza, I am newbie in openfoam and using reactingFoam in of5.
I am looking into solution for species source term using the same code you mentioned.
But, I have found the following error in attachment. Could you please tell me about any solution for this?

[Mirza8]

Hi Nuralam,

I do not have OF5 to test the same code right now. But the error seems to be related to the surfaceScalarField phi, which is missing!
I think you need to make sure that your modifications are made after the createFields.H, in which the

#include "compressibleCreatePhi.H"

is called to create phi.

Best wishes,
Morteza

[run]

Great. Now, its working. Thank you.

-NurAlam

[run]

Hi Morteza

I want to measure the following parameter:

dydt [n]= {Y(n-1) - Y (n+1)} / {t(n-1)-t(n+1)}

where, Y(n-1) and Y(n+1) are mass fraction of species at time t(n-1) and t(n+1) respectively.

Could you please tell me how I can store time date and calculate dydt?

Redards,
NurAlam