CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Post-Processing

Get the rate of reactions from the reactingFoam solver

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

Like Tree3Likes
  • 1 Post By Mirza8
  • 1 Post By Paulo SPS
  • 1 Post By Mirza8

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 5, 2019, 09:46
Default Get the rate of reactions from the reactingFoam solver
  #1
New Member
 
Morteza Mousavi
Join Date: Jul 2019
Location: Lund, Sweden
Posts: 4
Rep Power: 2
Mirza8 is on a distinguished road
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
Code:
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:

Code:
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 likes this.
Mirza8 is offline   Reply With Quote

Old   August 6, 2019, 10:10
Default
  #2
New Member
 
Join Date: Oct 2017
Posts: 10
Rep Power: 4
Paulo SPS is on a distinguished road
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):


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

Old   August 6, 2019, 11:19
Default
  #3
New Member
 
Morteza Mousavi
Join Date: Jul 2019
Location: Lund, Sweden
Posts: 4
Rep Power: 2
Mirza8 is on a distinguished road
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:

Code:
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?
Attached Images
File Type: png Screenshot_2019-08-06_15-24-52.png (138.0 KB, 4 views)
Mirza8 is offline   Reply With Quote

Old   August 7, 2019, 04:55
Default
  #4
New Member
 
Join Date: Oct 2017
Posts: 10
Rep Power: 4
Paulo SPS is on a distinguished road
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:


Code:
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 likes this.
Paulo SPS is offline   Reply With Quote

Old   August 7, 2019, 12:22
Default
  #5
New Member
 
Morteza Mousavi
Join Date: Jul 2019
Location: Lund, Sweden
Posts: 4
Rep Power: 2
Mirza8 is on a distinguished road
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:

Code:
// 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?
Attached Images
File Type: png source O2.png (55.1 KB, 2 views)
Mirza8 is offline   Reply With Quote

Old   August 13, 2019, 11:57
Lightbulb Possible solution!!!
  #6
New Member
 
Morteza Mousavi
Join Date: Jul 2019
Location: Lund, Sweden
Posts: 4
Rep Power: 2
Mirza8 is on a distinguished road
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 r_i, this code creates a scalar field which contains the production rate of species s_i through that reaction:

Code:
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, \sum_{r_i}RR(r_i, s_i), 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!
Paulo SPS likes this.
Mirza8 is offline   Reply With Quote

Reply

Tags
reactingfoam, reaction rate, source terms

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
How are chemical reactions stored after parsing? DanielW OpenFOAM Programming & Development 3 June 7, 2018 15:39
No reactions in reactingFoam 2.1 OMN OpenFOAM Running, Solving & CFD 16 April 7, 2015 13:14
mass fraction of the products in reactingFoam solver yash.aesi OpenFOAM Running, Solving & CFD 6 October 19, 2013 05:19
adding radiation to reactingFoam solver Marshak OpenFOAM 1 February 13, 2012 06:43
Finite rate chemical reactions Marc Segovia Siemens 3 November 28, 2001 05:21


All times are GMT -4. The time now is 06:02.