View Issue Details

IDProjectCategoryView StatusLast Update
0001348OpenFOAM[All Projects] Bugpublic2015-02-24 20:08
Reporteruser281Assigned Tohenry 
Status resolvedResolutionfixed 
Platformx86-64OSUbuntuOS Version12.04
Product Version 
Fixed in Version 
Summary0001348: Foam::constTransport<Thermo>::operator+= and -= crash when multiple species have 0 nMoles
DescriptionFoam::constTransport<Thermo>::operator+= and Foam::constTransport<Thermo>::operator-= will crash if the first two species of a reactingMixture has 0 moles. This is caused due to line 213 (+= operator) and 237 (-= operator) in constTransportI.H.
Steps To ReproduceCreate a foam case with a reactingMixture of 3 gases A, B, C. In the reactions list, set the order as {A, B, C}. In the initial setup of the problem, fill the domain with gas C and run the simulation. thermo.correct will loop over the gases in the order set in reactions and the program will crash when calculation Pr since both A and B will have 0 moles.
Additional InformationPotential fix for the problem in constTransportI.H.

Current line 213:
rPr_ = 1.0/(molr1/rPr_ + molr2/st.rPr_);

Suggested correction:
rPr_ = molr1/rPr_ + molr2/st.rPr_;
rPr_ = max(rPr_, SMALL)*sign(rPr_);
rPr_ = 1./rPr_

Current line 237:
rPr_ = 1.0/(molr1/rPr_ - molr2/st.rPr_);

Suggested correction:
rPr_ = molr1/rPr_ - molr2/st.rPr_;
rPr_ = max(rPr_, SMALL)*sign(rPr_);
rPr_ = 1./rPr_
TagsNo tags attached.



2015-02-07 17:45

updater   ~0003716

This bug report feels familiar... I did some back tracking of the changes to this file and:

 1. These operators were added to "thermophysicalModels/specie/transport/const/constTransportI.H" back on the date 2013-03-25, for fixing this bug report:

 2. Later on, date 2013-08-19, there was a correction for the addition of the Prandtl numbers, so that they would be added with the.. sort-of inverse values.

 3. And I know there is yet another bug report similar to this one on the need to invert the addition of the inverses... but I can't find it.

@LWhitson2: Any chance you can provide the test case you've mentioned? I ask this because I tried just now to set-up such a case based on the tutorial "combustion/reactingFoam/ras/counterFlowFlame2D" and I wasn't successful in reproducing the same exact error with OpenFOAM 2.3.x.


2015-02-09 13:19

reporter   ~0003734

Why would you want to set nMoles to 0?

This basically means that the component is not part of the mixture and can be removed entirety. Also saves you one transport equation to solve.



2015-02-09 15:47


Last edited: 2015-02-09 15:48

View 2 revisions

The basis of this calculation is a premixed flame. For example, let's assume the case where C + M -> A + B + M as our base reaction. Now, let's assume that the reaction is occurring in a co-flow medium (D). On the first time step, the code will loop through the species in the list in the order which they are specified. Since A and B both have 0 moles at the first time step (since they are both product gases that will only be prevalent after the initial combustion) the code will crash on the calculations of the various values. By filling the domain with gas C to start with, I am simplifying the problem a little bit with less variables just to show the problem. It would work the same way if you were to have 4 gases (A, B, C, D) and A, B were not prevalent at the start.


2015-02-09 16:21

reporter   ~0003740

I understand the implications for combustion simulation, but you are setting up your case in the wrong way!
By setting nMoles to 0 for species A, this species would never be taken into account in the thermodynamic property evaluation, no matter which value Y_A has during the simulation.

Just set nMoles to 1 for "single molecule" species, i.e. N2, O2, CO2, H2O, OH, ...
Then, in your 0 directory you initialise only species C and M with appropriate values (Y_C + Y_M = 1) and set Ydefault to 0. The other species mass fraction (here A and B) will be set to 0 for the initial time step.
Have a look at the tutorials that use multiComponentMixture or reactingMixture.


2015-02-09 19:15

manager   ~0003744

In the limit of having 0 moles of either species I think in the case of + it would be more logical to set rPr_ to that of one of the species rather than 1/SMALL. However in the case molr1/rPr_ - molr2/st.rPr_ = 0 it is not so obvious what the logical result should be.


2015-02-09 21:28

manager   ~0003745

I have push a fix to constTransport in OpenFOAM-2.3.x which should resolve the problem with only a very small overhead. Could you please test and let me know if there are any issues. Thanks.


2015-02-24 18:42


Henry, your code fixed the + and - operators, but the same error is still present in the += and -= operators. I believe you could fix those the same way, couldn't you?


2015-02-24 19:01

manager   ~0003898

I will check. I have also been considering the current method of averaging Pr rather than rPr; it would be much simpler to average rPr and I don't see any reason why we shouldn't and there is no physical basis for either approach.


2015-02-24 19:03


In the end I feel that it is a moot point because the end result should always be the value of the first non-zero species. That is why I used SMALL in my example because it was quick and easy.


2015-02-24 19:53

manager   ~0003900

In the case of several species each with non-zero nMoles how should the average Pr be calculated? Perhaps we should store the conductivity instead and average that buth then how should the conductivity or the viscosity be averaged? The current approach is very approximate.


2015-02-24 20:07

manager   ~0003901

Resolved by commit 59603afc8f1ef2427bd468820f6fcd229be9e6c1

Issue History

Date Modified Username Field Change
2014-07-11 14:44 user281 New Issue
2015-02-07 17:45 wyldckat Note Added: 0003716
2015-02-09 13:19 dkxls Note Added: 0003734
2015-02-09 15:47 user281 Note Added: 0003739
2015-02-09 15:48 user281 Note Edited: 0003739 View Revisions
2015-02-09 16:21 dkxls Note Added: 0003740
2015-02-09 19:15 henry Note Added: 0003744
2015-02-09 21:28 henry Note Added: 0003745
2015-02-13 12:33 henry Status new => resolved
2015-02-13 12:33 henry Resolution open => fixed
2015-02-13 12:33 henry Assigned To => henry
2015-02-24 18:42 user281 Note Added: 0003897
2015-02-24 18:42 user281 Status resolved => feedback
2015-02-24 18:42 user281 Resolution fixed => reopened
2015-02-24 19:01 henry Note Added: 0003898
2015-02-24 19:03 user281 Note Added: 0003899
2015-02-24 19:03 user281 Status feedback => assigned
2015-02-24 19:53 henry Note Added: 0003900
2015-02-24 20:07 henry Note Added: 0003901
2015-02-24 20:07 henry Status assigned => resolved
2015-02-24 20:07 henry Resolution reopened => fixed