
[Sponsors] 
Solver Validation chemFoam using analytical Solution 

LinkBack  Thread Tools  Search this Thread  Display Modes 
October 31, 2012, 06:35 
Solver Validation chemFoam using analytical Solution

#1 
Member
Join Date: Nov 2010
Location: Tokyo / Japan
Posts: 40
Rep Power: 12 
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 (ACell, also open). To understand chemFoam, I implemented the reaction Code:
A + B > C 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] the Molecular Masses to be Code:
M_A = 1 [g/mole], M_B = 1 [g/mole] Code:
M_C = 2 [g/mole] ACell 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 Allrunscript. 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 Last edited by Hanzo; October 31, 2012 at 23:06. 

November 5, 2012, 21:51 

#2 
Member
Join Date: Nov 2010
Location: Tokyo / Japan
Posts: 40
Rep Power: 12 
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 

November 6, 2012, 04:20 

#3 
Member
Michiel
Join Date: Oct 2010
Location: Delft, Netherlands
Posts: 97
Rep Power: 12 
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?


November 7, 2012, 04:01 

#4 
Member
Join Date: Nov 2010
Location: Tokyo / Japan
Posts: 40
Rep Power: 12 
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 } 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 00098 virtual tmp<volScalarField> rho() const 00099 { 00100 return p_*psi(); 00101 } http://www.cfdonline.com/Forums/ope...openfoam.html ) Actually, in the program ACell 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 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 ACell? Last edited by Hanzo; November 7, 2012 at 04:22. 

October 26, 2018, 10:47 

#5 
Member

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.


April 4, 2019, 09:08 

#6 
New Member
Join Date: Mar 2019
Posts: 4
Rep Power: 3 
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 secondorder 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 ACell, however the fundamental process I use in Mathematica is the same, i.e. a nontemperature 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). 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
analytical solution of 1D shock tube problem  ma  Main CFD Forum  10  May 11, 2015 06:20 
thobois class engineTopoChangerMesh error  Peter_600  OpenFOAM  4  August 2, 2014 10: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 