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 Code:
constant k=3 [m^3 / (mole sec)], the Molecular Masses to be Code:
M_A = 1 [g/mole], M_B = 1 [g/mole] Code:
M_C = 2 [g/mole] 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)] I figuered: If I reduce k by a factor of 7: Code:
=> k =428.5 [m^3 / (kmole sec)] Code:
=> k =4.285e+5 [cm^3 / (mole sec)] 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 |
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 |
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?
|
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); 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() basicPsiThermo.H: http://foam.sourceforge.net/docs/cpp/a07886_source.html Code:
00097 //- Density [kg/m^3] - uses current value of pressure 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) 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? |
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.
|
2 Attachment(s)
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). |
All times are GMT -4. The time now is 04:11. |