[Hanzo]

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

[Hanzo]

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

[michielm]

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?

[Hanzo]

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:

ODEChemistryModel.C
http://foam.sourceforge.net/docs/cpp/a07939_source.html

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

thermo()->rho()

and this is computed using:
basicPsiThermo.H:
http://foam.sourceforge.net/docs/cpp/a07886_source.html

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
http://www.cfd-online.com/Forums/ope...-openfoam.html
)

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?

[mkhm]

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.

[rando_foamer]

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).