CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Solver Validation chemFoam using analytical Solution (http://www.cfd-online.com/Forums/openfoam-solving/108728-solver-validation-chemfoam-using-analytical-solution.html)

Hanzo October 31, 2012 06:35

Solver Validation chemFoam using analytical Solution
 
3 Attachment(s)
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

Code:

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

Code:

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
Code:

M_A = 1 [g/mole], M_B = 1 [g/mole]
and therefore
Code:

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
Code:

k = 3 [m^3 / (mole sec)] => k =3000 [m^3 / (kmole sec)]
and using the chemkinReader
Code:

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:
Code:

=> k =428.5 [m^3 / (kmole sec)]
or for the chemkin Case
Code:

=> 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
http://www.cfd-online.com/Forums/att...1&d=1351679335

Fig 2: k / 7

http://www.cfd-online.com/Forums/att...1&d=1351679326

Hanzo November 5, 2012 21:51

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 :o

michielm November 6, 2012 04:20

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 November 7, 2012 04:01

1 Attachment(s)
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
Code:

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
Code:

thermo()->rho()
and this is computed using:
basicPsiThermo.H:
http://foam.sourceforge.net/docs/cpp/a07886_source.html
Code:

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

Code:

    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 :confused:

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

http://www.cfd-online.com/Forums/att...1&d=1352275184

maybe I should somehow use the density given by openFoam to convert the results given by analytical solution and A-Cell?


All times are GMT -4. The time now is 05:30.