
[Sponsors] 
October 31, 2012, 06:35 
Solver Validation chemFoam using analytical Solution

#1 
Member
Join Date: Nov 2010
Location: Tokyo / Japan
Posts: 40
Rep Power: 6 
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: 6 
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: 6 
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: 6 
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 } How rho_ and psi_ are calculated in compressible solvers of OpenFOAM? ) 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. 

Thread Tools  
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 05:20 
thobois class engineTopoChangerMesh error  Peter_600  OpenFOAM  4  August 2, 2014 09: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 