CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Solver Validation chemFoam using analytical Solution

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

LinkBack Thread Tools Search this Thread Display Modes
Old   October 31, 2012, 06:35
Default Solver Validation chemFoam using analytical Solution
Join Date: Nov 2010
Location: Tokyo / Japan
Posts: 40
Rep Power: 12
Hanzo is on a distinguished road
Hello Foamers,

the last days I investigated the chemistry functionalities of OpenFoam 2.0.x and I created a simple test case for chemFoam, computed the analytical solution and compared to another chemistry solver (A-Cell, also open).

To understand chemFoam, I implemented the reaction

A + B -> C
to compare with the analytical solution I chose the following conditions:

constant k=3 [m^3 / (mole sec)],
A0= 0.5 [mole/m^3], B0= 0.3 [mole/m^3] and C0= 0 [mole/m^3],
0 [sec] <= t <= 10 [sec]
For conversion from Mass fraction back to moles after using chemFoam, I consider
the Molecular Masses to be
 M_A = 1 [g/mole], M_B = 1 [g/mole]
and therefore
M_C = 2 [g/mole]
The Situation I have:
A-Cell and the analytical solution coincide perfectly. The Solution given by chemFoam reaches the steady state way too fast (Figure 1).

I computed results using the native foam chemistry reader using
k = 3 [m^3 / (mole sec)] => k =3000 [m^3 / (kmole sec)]
and using the chemkinReader
k = 3 [m^3 / (mole sec)] => k =3e+6 [cm^3 / (mole sec)]
Good news about that: the results differ exactly the same from the analytical solution.

I figuered: If I reduce k by a factor of 7:
=> k =428.5 [m^3 / (kmole sec)]
or for the chemkin Case
=> k =4.285e+5 [cm^3 / (mole sec)]
The error is of the order where I could start to argue that it depends on the ode solver I am using (just tried ODE solver SIBS) (Figure 2).

Why do I have to reduce the rate constant k by this strange factor around 7 to achieve meaningful results? Please have a look on the figures I append.
Any ideas or comments are appreciated.

Also, I append my case, with an Allrun-script. Please have a look, it should run fully automatic and one can play easily with the value k in chemkin/chem.inp or constant/reactions

Fig 1: regular converted k

Fig 2: k / 7

Attached Images
File Type: png k_regular.png (6.9 KB, 216 views)
File Type: png k_by_seven.png (6.9 KB, 217 views)
Attached Files
File Type: gz simple_reaction.tar.gz (70.3 KB, 21 views)

Last edited by Hanzo; October 31, 2012 at 23:06.
Hanzo is offline   Reply With Quote

Old   November 5, 2012, 21:51
Join Date: Nov 2010
Location: Tokyo / Japan
Posts: 40
Rep Power: 12
Hanzo is on a distinguished road
Sorry for pushing this up, but I still have no new clues on my problem

Any suggestions how I can explain my problem better or any ideas are welcome
Hanzo is offline   Reply With Quote

Old   November 6, 2012, 04:20
Join Date: Oct 2010
Location: Delft, Netherlands
Posts: 97
Rep Power: 12
michielm is on a distinguished road
This sounds a lot like a unit conversion issue. I guess you've already done this, but just to be sure: have you checked that all units are consistent?
michielm is offline   Reply With Quote

Old   November 7, 2012, 04:01
Join Date: Nov 2010
Location: Tokyo / Japan
Posts: 40
Rep Power: 12
Hanzo is on a distinguished road
Hi Michiel,

thank you very much for your answer.

I rechecked the Units but still I think they are correct ...

However, I would like to give an update on my problem:

I looked inside ODEChemistrySolver.C and I saw that the concentration Ci are calculated using the following:

00721             scalarField c(nSpecie_, 0.0);
00722             for (label i=0; i<nSpecie_; i++)
00723             { 
00724                 const scalar Yi = Y_[i][celli];
00725                 c[i] = rhoi*Yi/specieThermo_[i].W();
00726             }
The resulting c[i] are then of [kmol/m^3 ] and these are plugged into the selected ODE Solver.

I was a little bit confused how the density rhoi comes into play and I figured its and inherited property of the class
basicPsiThermo. Then rho is given by
and this is computed using:
00097             //- Density [kg/m^3] - uses current value of pressure
00098             virtual tmp<volScalarField> rho() const
00099             { 
00100                 return p_*psi(); 
00101             }
(Other users also discussed about that

Actually, in the program A-Cell and in the analytical solution there is no density defined.
c[i] just comes from the user and is given in [mole/m^3] but in OpenFoam density plays a major role.

The final value of rho in my simulation is 0.0075 kg/m^3 which is 7.5 g/m^3 and, voilą, I found a factor similar to the value of 7 I mentioned in my first post.

So I adapted ODEChemistryModel.C defining

     kf = R.kf(...) / (1000*rho)
The result looks much better than in the beginning (see figure 3) but I am totally unhappy about this
solution because it is not physical (units make no sense), and there is still a
difference to analytical solution

figure 3: k divided by density given in [g/m^3]

maybe I should somehow use the density given by openFoam to convert the results given by analytical solution and A-Cell?
Attached Images
File Type: png update01_k_by_1000_times_rho.png (7.1 KB, 198 views)

Last edited by Hanzo; November 7, 2012 at 04:22.
Hanzo is offline   Reply With Quote

Old   October 26, 2018, 10:47
Join Date: Jul 2017
Posts: 87
Blog Entries: 1
Rep Power: 5
mkhm is on a distinguished road
Did you get any answer to your question ? I have the similar problem. By looking at ChemistryModel, it seems that the unit of reaction rate coefficient should be m3/kmol/s but when I consider this unit for my reaction rate coefficients, I have inconsistent results compared to the analytical one.
mkhm is offline   Reply With Quote

Old   April 4, 2019, 09:08
New Member
Join Date: Mar 2019
Posts: 4
Rep Power: 3
rando_foamer is on a distinguished road
I realise this thread is ooooold and the original poster has probably moved on to bigger and better things, but I recently had the same problem and managed to sort out the issue. Since this thread helped me massively with understanding how chemFoam goes about things, I thought I'd post my solution so that others with the same problem can get past it.

Firstly, chemFoam assumes that the units for rate constants provided in the chemkin chem.inp file are in cm^3 mol^-1 s^-1 for second-order rate constants. In other words, all constants should be provided with units of cm for length, mol for quantity, and s for time.

It carries out its own conversions of the rate constant, but we don't need to concern ourselves provided the other parameters of the system are correct. In the initialConditions file, the pressure supplied should be appropriate to the number of moles in the system and the temperature.

In Hanzo's case setup, the pressure he/she specifies (12000 Pa) does not match the number of moles (0.8 mol) and the temperature (310 K). By my calculation, it should be ~2062 Pa.

The pressure is the only way that chemFoam can work out the actual concentrations in the mixture, since whether fractionBasis is given as moles or as mass, chemFoam immediately converts these numbers to mass fractions. So, if the pressure is not correct, the density is not correct and the concentrations calculated during simulation will not be correct. This is why Hanzo's original chemFoam simulated rates that were much too high - the pressure provided suggested to chemFoam that the concentrations were much higher than intended.

Side note: I would also suggest that the constantProperty in the initialConditions file should be set to volume rather than pressure, as in this case the pressure will change throughout reaction as the number of moles does not remain constant, whereas the volume does remain constant (1 m^3).

I ran this through chemFoam and compared it against equivalent simulations in Mathematica (unlike Hanzo, I don't use A-Cell, however the fundamental process I use in Mathematica is the same, i.e. a non-temperature aware set of ODEs). The two now agree - see the figure.

I have also attached the case setup I used for the chemFoam simulations (using OpenFoam v1812). It can be run using the "run" script, and will output probes for AL, BL, CL, p, and T in the postProcessing folder. Note that AL, BL, and CL will be given as mass fractions, not concentrations! You have to convert back to moles or concentrations (or equivalently, convert the results from other simulations to mass fractions).
Attached Images
File Type: png comp.png (11.4 KB, 6 views)
Attached Files
File Type: gz Hanzo.tar.gz (30.0 KB, 2 views)
rando_foamer is offline   Reply With Quote


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
analytical solution of 1D shock tube problem ma Main CFD Forum 10 May 11, 2015 06:20
thobois class engineTopoChangerMesh error Peter_600 OpenFOAM 4 August 2, 2014 10:52
Analytical solution to Euler equations Hooman Main CFD Forum 1 November 5, 2012 04:39
Logging files (.C and .H) the solver uses during the solution of a case Bufacchi OpenFOAM 2 February 1, 2012 09:07
Validation of axisymmetric Euler solver Fernando Velasco Hurtado Main CFD Forum 3 December 19, 1999 19:59

All times are GMT -4. The time now is 03:15.