I tried to introduce kinetic dependence into Saturne 7.0.2 model but without any positive results. As I understand, from three gas combustion models available, only Libby-Williams one deals with kinetics, but it does not work correctly. I investigated the code in extent I can and found the following.
Main part of Libby-Williams model in Saturne is in
pdflwc.f90 file. In this file component mass fraction, mixture properties (including temperature) and kinetic reaction rate (form of Arrhenius law) is calculated with or without PDF (statistical distribution function) depending on concentrations statistical dispersions (to save calculation time if dispersion is very low). Common algorithm in this file,
applied to each mesh cell, is:
1. Initialize, set pointers...
2. Get mixture enthalpy
h in cell:
Code: Select all
h = ((hmax-hmin)*f + hmin*fmax - hmax*fmin)/(fmax-fmin)
where
f is a sort of reaction progress variable or reacted fuel fraction (sorry, the code is very difficult to understand, most variables are not commented including
doxygen and has not "clear" names like
fuelMassPart or so, they named like "f", "y" etc, it's not easy for the user to get what they mean from the code).
I suppressed array indexes to make expressions more suitable to get their meaning, the problem is not in PDF / Dirac's function etc, I tried without PDF.
3. Calculate components (fuel, oxydant and products) mass fractions (
yfuel etc) in cell:
Code: Select all
yfuel = y
yoxyd = 1 - (coeff3+1)*f + coeff3*y
yprod = 1 - yfuel - yoxyd
yo2 = coeff1 - (coeff1 + coeff2) * f + coeff2 * y
where
coeff3 and coeff1 are stoichiometric coefficients or something like this, the problem is not in stoichiometry.
4. Get mixture temperature and density in cell by enthalpy (
cothht routine with mode=1).
5. Calculate kinetic reaction rate
w in terms of the fuel consumption per second. It's a form of standard Arrhenius law + multiplication of fuel and oxygen (not oxydant) mass fractions, it's very common and OK, the problem is not in kinetic rate itself.
6. Set mass fractions, temperature, density and kinetic reaction rate (fuel mass source) in program-wide arrays.
7. Free memory/pointers etc.
Other important part of the model is in
lwctss.f90 file. In this file, fuel concentration field sources are calculated. It's very important that
the fuel field is associated with kinetic rate w, but other component fields are not.
The problem, in my understanding, is that product concentration (
yprod), enthalpy and temperature are governed not by the reaction rate but by reaction progress variable and, maybe fuel concentration (
y and
f variables, sorry, it's not clear what they mean). If you disable
lwctss subroutine (place return in it's beginning), you will see that reaction progress variable and fuel concentration are solved as transported scalars, but oxydant and product mass fractions are not, they are just a state functions in mesh cells. So even if initial/inlet temperature is very low (30*C), model fuel reacts fine! You cannot affect it by any kinetic parameters in input data. In the same cell (in Paraview), reaction rate (chemical source) tends to zero (10^-13), but temperature increases. The model behaves like it's high kinetic rate while it's zero (anyone can see it in my complete example in first post, just run it). As we know from reaction basics, if diffusion or reaction rate is zero, summary rate must also be zero, but in calculation, in opposite, it's very high.
The mechanism of this, as I can see, is follows. The program calculates progress variable and fuel concentration like transported scalars. Due to turbulent dispersion/mixing we have some intermediate (not 0 or 1) values of these scalars in cells. Then it calculates oxydant and fuel concentrations that are, roughly, arbitrary due to turbulent dispersion (or some balance error arises from reactant mixing/turbulence). Products concentration is calculated by balance:
yprod=1-yfuel-yoxyd. So we have some product concentration (example in cell near inlet:
ymfuel->0, ymoxyd=0.894, ymprod=0.106, T=447K, w->0; at inlet:
Tin=300 K). Then, with reagents mixing downstream, product concentration increases up to unity in some cells and temperature gets higher up to adiabatic (in simple tests I found that enthalpy and temperature depend on
yprod). Thus, we have fast reaction with kinetic-diffusion model at 30*C initial temperature that is non-physical. Kinetic part only affects fuel field source and cannot stop this. Even if you set
w(idirac)=0 directly in the code you change quite few in test results.
I tried to correct this by setting yprod to zero at low temperatures etc but without luck. There are some artifacts and no correct behavior. The problem is that model based not on complete diffusion-kinetic rate and corresponding component sources like in standard approach but on implicit variables like "progress variable" not governed
directly by the real reaction rate. So I didn't find an option to fix this. From my point of view, complete new kinetic-diffusion gas combustion model is required that operates with kinetc-diffusion rate of a given set of reactions and real component fields like CH4 and O2, not progress variable or products mass fractions. Such approach works well in CFX and Fluent, we also used kinetic-diffusion model for coal (not gas) combustion in our engineering programs for decades. An astronomic level is CHEMKIN import and complex solver for detailed mechanisms like GRI-Mech-3.0, but it's extremely hard to implement (Fluent can do it but with problems, not for every task). That's why an optimal solution is simple diffusion-kinetic algorithm with real components fields (CH4, CO, O2, CO2, N2, H2O). Maybe you will introduce it in future? For now, because I can be wrong, please, approve that current model cant take into account kinetic (flammability dependence from temperature) and, maybe, there is a workaround that I didn't find? I cannot introduce complete new combustion model in Saturne, it's outside of my knowledge... I can implement some engineering method described by algebraic relations but not complex things like numerical turbulent diffusion solving, PDFs, flamelets or so...