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 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: 9
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, 130 views)
File Type: png k_by_seven.png (6.9 KB, 131 views)
Attached Files
File Type: gz simple_reaction.tar.gz (70.3 KB, 15 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: 9
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: 9
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: 9
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, 114 views)

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


Thread Tools
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 On
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 05:20
thobois class engineTopoChangerMesh error Peter_600 OpenFOAM 4 August 2, 2014 09: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 12:14.