Skip to content

GitLab

  • Menu
Projects Groups Snippets
    • Loading...
  • Help
    • Help
    • Support
    • Community forum
    • Submit feedback
    • Contribute to GitLab
  • Sign in / Register
  • openfoam openfoam
  • Project information
    • Project information
    • Activity
    • Labels
    • Planning hierarchy
    • Members
  • Repository
    • Repository
    • Files
    • Commits
    • Branches
    • Tags
    • Contributors
    • Graph
    • Compare
  • Issues 379
    • Issues 379
    • List
    • Boards
    • Service Desk
    • Milestones
  • Merge requests 13
    • Merge requests 13
  • Deployments
    • Deployments
    • Releases
  • Wiki
    • Wiki
  • Activity
  • Graph
  • Create a new issue
  • Commits
  • Issue Boards
Collapse sidebar
  • Development
  • openfoamopenfoam
  • Issues
  • #1402

Closed
Open
Created Aug 19, 2019 by Admin@OpenFOAM-adminMaintainer

Possible mistake in omega calculation in pyrolysisChemistryModel.C

Hi,

I was trying to use solidChemistry and pyrolysisChemistryModel in a solver, and I was getting weird results or runtime errors (undefined behavior) for a simple case where reactions took place in a solid medium! When I tracked down the runtime error, I found that the problem arises from the line 267 of pyrolysisChemistryModel.C, where we have

const scalar exp = R.lhs()[si].exponent;

where we get strange values for the "exp" variable! If you take a look at the whole block of this code:

const label Nl = R.lhs().size();

for (label s=0; s<Nl; s++)
{
    label si = R.lhs()[s].index;
    const scalar exp = R.lhs()[si].exponent; // This is line 267, the problematic line IMO

    kf *=
        pow(c1[si]/Ys0_[si][celli], exp)
       *(Ys0_[si][celli]);
}

return kf;

I believe the "si" at line 267, should be replaced by "s", so the new code should be:

const scalar exp = R.lhs()[s].exponent;

I tried this and it solved the problem in my code! Because as far as I understand, "s" is the species index in the R.lhs() (left hand side of the reaction) and "si" is the species index in the species list! In case it is not clear or it is not convincing, please take a look at the following reaction file as an example:

species
(
	A
	B
	C
);

gaseousSpecies
(
	D
	E
	F
);

reactions
{
	reaction1
	{
		type        irreversibleArrheniusSolidReaction;
		reaction    "B = C";
		A           1.5e14;
		beta		0;
		Ta          23653.7493709109;
		Tcrit       500;
	}
}

So in this case, to calculate the omega of the reaction, "s" will be 0 but "si" will be 1! Therefore, trying to read R.lhs()[si].exponent leads to undefined behavior, because R.lhs() has a size of 1!

I tried to explain the problem as clear as I could, so sorry for the long post, and please let me know if there is still something missing.

Best regards, Morteza

Assignee
Assign to
Time tracking