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)]`
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,

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.

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.